{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1999~2005的预测值： [2.67       3.11608066 3.25719252 3.40469464 3.5588764  3.72004029\n",
      " 3.8885025 ]\n",
      "\n",
      "-------------------------------\n",
      " 相对误差 [ 0.         14.42307692 17.84615385 21.70087977 25.         28.22580645]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyvUlEQVR4nO3dd3hU1brA4d+aJJOEVEgCoYcSEJQiRIrSAoIoiCBcz5FjQe8BFVGvjWs/etRjARtNmiAooIg0jYKIgPQmzQIiEGogIQkhdZLMrPvHDLlJSJlAkp2Z+d7nmSdklb2/lQlf9qy99t5Ka40QQgjPYDI6ACGEENVHkr4QQngQSfpCCOFBJOkLIYQHkaQvhBAexNvoAEoTHh6uo6KijA5DCCFcyu7du89rrSNKq6+xST8qKopdu3YZHYYQQrgUpdTxsuplekcIITyIJH0hhPAgkvSFEMKDSNIXQggPIklfCCE8SI1dvSOEEB5lQjRkJl5eHlAXnj1cabtx6aRvsVhISUkhPT0dq9VqdDiiELPZTHh4OCEhIUaHIoRrKCnhl1V+hVw26VssFk6cOEHt2rWJiorCx8cHpZTRYQlAa012djanTp3C19cXPz8/o0MSQji47Jx+SkoKtWvXJjw8HLPZLAm/BlFKUatWLcLDw0lKSjI6HCFqvFxLTrXty2WTfnp6OsHBwUaHIcoQFBRETk71/TIL4Woslmy2LZ5IylvXVts+XXZ6x2q14uPjY3QYogze3t7k5+cbHYYQNU++hV9WTKH+gWl04zx/+rQmMu98tezaZY/0AZnSqeHk/RGiqJzsLGw7PoFJneh04N+k+4Tza+xcop/fZl+lU5LSyq+Qyx7pCyGEq8jOyuKXlZNpeXAG9UiGRl2wDv6I6JZ9USbHsXclLsssiyR9F5CRkYGfnx/e3vJ2CeFKsrIy+WX5ZFr+OZObSOaQT1vS+35Iy26342XQJ2GXnt7xBPHx8YSHh7N06dJq3e/Zs2cvK0tOTiY3N7da4xDCJeXlwI5ZZE1sT48/3+KiuR4H+39G6xe20LL7EDBw6lOSfg0XFRVF+/btmTlzZrlt+/Xrh1Kq4OXl5cW5c+fo06dPkfLCr7179162nQ0bNtCwYUPWrl1Lfn4+WuuC7Y8ePbqgnc1mIy8vr9LGKoSry8jMYNPC/2Cb1BG+ewafsKb8ectntHp+C9fcZGyyv0SSfg2ycOHCEhPzzp07Wbt2bYl1CxcuLOhvNpt5+OGHSU1N5emnn6Zbt27Uq1ePVatWYbVa0VoXeVksFtq3b39ZHFOmTOHmm2/mwIED+Pj4YDKZUEqxb98+5s+fX+SPSuE/AkJ4qvSMdDZ89gZZE66jx5/vcMHcAO5bQcjYtbQy+Mi+OEn6xSRezOGuGVtJTK/+9eW+vr4AlyVnrTUZGRlYLJaC71NTUwGKXO1qNpvx9fUlNDSUNWvWMHTo0II2JtPlb7XZbL6sfOPGjSxZsoRXXnmFcePGkZeXh81mQ2tNnTp1+OSTTwpiyM3NdeoTiBDuypqbzfr5r5M9sR29j0wg1bcRf926kDrj1kLzPjUq2V8iSb+YSWsPszM+hUk/Vs+Z9MK8vLxo2rQpANOmTeM///kPFosFgLVr1xISEsLhw/a4vL296d27d5EL1C4tkTx69Cj79+9n2LBhJCYmcv78eS5cuHDZKzk5mdOnTxfctygzM5MxY8YA0LBhQ7y9vZk2bRrdunWjR48eZGVl8f7779OjRw+6dOnCihUrMJvN1fbzEaKmyM3OhG3T8ZrUkT5HJ5Lq15ijt31B6+c20rLroBqZ7C9xy+Ugf5ux9bKywe3rc2/3KLJzrYyau+Oy+t3HU8m36YLvP99+gs+3n0Ap6BJVh3u6NeX2Dg04cyGbJ7/ce1n/0T2bc3PbehxJyqBFROAVxT106NCCo/Po6GjGjx/P7NmzmTRpEiaTCYvFQpMmTQAIDAxk/fr1JW5nyZIlXHfddbRs2ZKOHTuyb9++Mvd77NgxmjRpwqhRo8jKyipSd/78eerXr8/y5cuLlPfp04eLFy9e0TiFcFWpF9LYs/xD2sXPJYJUaHoTuUNn0rplb6NDc5oc6Ts8f+s1DOnYAD9v+4/EpCAswEzHRqHVGseMGTN44IEHaNu2Lbt27eLhhx8mNzeXrKwsIiMjC6aAytK2bVtOnjxJQkICW7duLTgZe+mE7Lp169BaY7VaycjIoGnTpsTHx/Pzzz8ze/bsItsymUzExcURGhpa5LVp06YSp4yEcEcpF9JYO/dV8j9oT9/490n1b0LS8K/hge8wu1DCBzc90v/yoe6l1vmbvUqtf3HZASxWG77eJnKtNm69LpI3hrUr0qZBqH+Z27/So/xLevTowbfffkubNm2YNGkS48ePB+DDDz+kRYsWTm1j8ODBxMTE8MYbbzB16lSOHj1KWFhYkdscJyUlkZOTQ+PGjQFo3rw5f/zxR4mJfNCgQSUe6Qvh9vKyyd46G+vaifRTF/izVkeyb36BVp1vMTqyK+aWSf9Knc+w8I+uTRnZpQkLd5wgqZpP5mZmZnLNNdfwzTffMH/+fBo0aFBQFx8fT/PmzYu0z83NJS8vj4CAgMu2NWzYMN566y2mTp3KmDFj8PLyYvXq1QX1gwcPpmHDhkXW/9epU4cLFy5ctq1Vq1YRGRlZpCwlJYVRo0Zd4UiFqNmSUlI58+M0Ohz/FP/MRNLCbuBUz/G0un6A0aFdNUn6hcy4N6bg328Mva7a9x8YWP6nhPnz5xf5/tLUzCV5eXkcOHCAn376ieDgYA4dOsRPP/3E2rVri/R74YUXGDp0KNu3b6dr165l7nPgwIFypC88QmJKCnuXfsD1J+fRQaVhadwD3//6lMiom4wOrdJI0q9BTp06hZ+fH15eXkXK586dy3/+8x+2bdtGcHAwPj4+BUsmbTZbkbZKKW677TZ8fHxYtGgRb7zxBn369CE2NrZIuzvuuIPu3bvz3HPPsW7dujLjujSnX1hGRoYc6Qu3kZKayq6vJ9Lp5GcMUGn8GdiZvAEv0qBDP6NDq3RyJq4GadiwIWFhYQUnS4OCgpg7dy4vvvgic+fOZfLkyXTv3p1Vq1YRGhpKvXr1qF+/fkH//Px8vL292bZtG4cPH8bHx4cFCxbw0ksvlbi/F154gfXr1/Pjjz+WGdegQYMuW+7Zo0ePSh27EEaw5WTAlsmEzophwKkppARGk3DnMlo9+5NbJnyowJG+UmoocKfj26+01t9UZEdKKT/gV0ff5yvS19PEx8ezePFiZs+eTXp6OsuWLeOWW26hW7duBAYGMnr0aN59913eeecd+vfvX9Dv0i0RGjZsCNiPxseMGUPfvn2Ji4srWLp5aW39oEGDWLRoET179rxsG2lpaZw9e5acnBzy8vIum+vPz88nMzOTM2fOEBYW5tSqIiFqijNJ59m37D1uPLuAEFsapuZ9yOr+DK2ie5bf2dWVdPVn8RcwFkgFJgCzABtwnzN9C23jX0ACEORM+86dO+uy/P7772XWu6KlS5fqqKgoDejmzZvrt99+W1+4cOGydqdPn9Z33nmnBvT7779fUB4bG6ufeOKJErc9adIkHRUVpe+++26dmZlZagwnT57UgJ48ebL29fXVQUFBOiQkpMRXUFCQ9vPz03v37i11e+74PgnXdfJsko77+Dl9/pVGWv8rWP85oZ/OObLJ6LAqFbBLl5Fblda6hD8F/08pFQScBoZrrdc4ymYCbbTWTv1ZVEo1A34HxmmtP3GmT0xMjN61a1ep9X/88Qdt2rRxZlMuQ2vNlClTuOmmm+jUqVO57desWUOPHj3w9/evhuiujDu+T8IFWTKIX/0Rwbs/po5K53BQF0JvfZmItr2MjqzSKaV2a61jSqt3ZnrHF3uyXlOoLAHoWIE4PgIOAnMr0MfjKKV47LHHnG5feGpHCHG5EwmJ6B2zaHroE6Kykjka2g3bwJeIbuMB0zilKDfpa63PAwXrBJVS4cA9wBxndqCUGgjcDvwIfK6U+gN4X2udeUURCyFEOeLPJPLr8onceG4BdVQGuuXNqN7P0bzxDUaHZrgKLdlUSk0D/gb8DLzrZLe3HF/rOfZ3FzBcKdVda51dbPtjgDFAwT1mhBDCWcfPnOPAsoncmLiQwSqDwyHdYdDL1GntPuvsr1ZF1+lvBloCsUA3YGNZjZVSN2CfBpqltR7jKOsF/ATcD0wv3F5rPROYCfY5/QrGJoTwUDrnImrnLOpvnETT3Av8FdodNegVolvdaHRoNU6Fkr7WegGwQCn1OTCZ8uf1ox1f3y+0jZ+VUoeBDhXZtxBCFHf4xBn+WDGR/mlL8M9Pwyd6AGldn6Zly25Gh1ZjlZv0lVJmoJ7W+mSh4jhghBPbz3B8PVasPAeQh60KIa7IoeOnObhiIr2SFzNEZRBfpwdRw/+NatiZkPK7ezRnrsiNBQ4opeoUKosGjjvRdxegKXRUr5SqC7QGdlYgTo+WlZVV4rNo8/PzC+5/f+zYsYIHrDhj5cqVfPzxxyXWyUPRRU1x2ZPsci6yY97z1JtzA3ekzCGlTkfS7/2BqMfjoGFnY4N1Ec5M7/wEnALilFKvAZHAM8DzAEqp1oBFax1fvKPW+oxSahHwmVLqOcACvAycA5ZUygjcSKdOnTh79myRp1G98sorvPnmmyQkJGA2m0lLSyMwMBAvLy/y8/OJjIzkr7/+YubMmUyePJmPP/6YmJgYnn/+8ouev/zyy4IrZ3fu3Mm+fft45JFHirTZsGEDffv25YcffqB37954eXmhlKJfv3506NCBefPmAfaHolutVnx8fKrwJyI81oRoyEykLrAY4L1LFYouaI7U6Yn34Fdo0aKLYSG6qnKP9LXWecAgIAn7z/9l4AWt9VRHkxnAq2Vs4kFgBTAJ+Ar7tM6tWuvqfwhtDRcXF8eePXvYs2cPd999N/7+/txyyy0cOXKErKwsTp8+DcAvv/zChQsXyMjI4K+//gLgrbfe4vXXX2fKlCmkpKSwZcsWRo0axahRo7jttttYsWIFPj4+vPTSS/z+++94eXmVeO98eSi6qBEyE0up0DBmPS2e+JZASfhXxKkTuVrr48CQUur6lNPXAox3vEQZ6tevzxNPPEHTpk2ZM2cO27dvL7iPDtincAICAkp9mMqTTz7J448/zoEDB6hVqxa9e/dm48aNBTdHM5lMzJs3j759+5bY/9JD0Tdt2kTXrl0ZN25cwZF+WFgYEyZM4MEHHwTs9+gp72puIapEg+uNjsClyV02L5kQDa+GXP6aEF1+30rUu3dvXnzxRe6//36ioqIA+PbbbzGZTPTo0QNvb2/q1KlDaGgoJpOJb7/9Fq0177zzDqmpqUVuy5ySksIdd9zB+fPnC8pKO8KXh6ILo+VbbegLJ+H754wOxa3J/fQvKe3jZKkfMyuX1pqcnByGDRvGypUrOXPmTEGdr68vTZo0IT4+nqNHj6K1pkWLFkRFReHn58fFixdZsWIFM2fOJC4urqDfpXvyFL8XfnE2m00eii4Mk5Wbz/frN2HePplBegMK+QRZldwv6X//HJw9ULnbnDuoYu0j28Gtb1eoy4kTJ4iOjsbf35/MzEzMZjOPPvooeXl5rFixoqDdhAkTSElJ4csvvwTA29ubkJAQ1q9fz4gRI/j555/p0sU+16mUAij3hmyXHor++eefM2DA/z8OrvBD0QuTB6iIypCSmUvcmh+ou3cqQ/VW8pUP59uMpO4tz8CH7crfgLgi7pf0XVTTpk0LlkT26dOHUaNGYbVamTt3Ln5+fgXtsrOzS5zTN5vNLF++HJPJxJ49ewC4ePEiZrO53BU28lB0Ud1sx7fxx7wXuNe2m2xVi8R2D1P/lqeoG1jX3iCgbsmfsgPqVm+gbsj9kn4Fj7ALvFrGJR0PxJVeV4UiIiIYN25ckbK0tDRiYkq+a2pSUhJjx47loYceAiAhIaHIw9UtFkup+5KHoouq9vvpNHav/5p7cpdgOr6ZG3xrc77DeMJjH8XfP7Ro42edv+ZEVIz7JX03MmSIfcFU4YeaHzx4sOCEK1Dkoq1Zs2axfft2XnrpJQYMGMD27dtp3LgxgYGB3H///WitsVqtJV7oVRp5KLq4Glprth5JYteqz+iT+Bn3mo6RFxCJz8C3MXe6j3BzgNEhehxJ+pfU4I+TOTn2SxqOHDnC4cOH6dzZfuVho0aNClbrWK1WZsyYwTPPPMP111/PzJkz6dWrF7GxsezYsYNPP/0UgNzc3DKP+IuTOX1xpc6lprPwkw8YfPELHjed5kJAY7J6f0itmJHgLY/XNIok/Utq8MfJQYMGcdtttzF8+HAGDhzI9u3bOX78OJs2bSpo88knn2CxWAo+BWzdupWNGzcyZswY7rjjDubMmcO9997LuXPnACjrqWTF9y1H+sJZOXlWjpw5z7WJ31B380c8mXGC1JBocvvNIrTdneAlKcdo8g7UIAkJCZw/f55Tp05Rq1atgvL09HQeffRRdu7cybZt20hOTuaVV15h9erVzJs3j7y8PF566SXGjh1LrVq1yM3N5amnnmLEiBHcc889JCcnM3bsWHr06EGzZs0ASpzikYeiiyuVlp3H4s2/k7VlFv+wrQTSUI1ugFsnULvVLeBYSSZqgLIeoGvkyxMfjL5y5UodERGhBw4cqFNSUrTWWn/xxRc6MjJSDxgwQJ89e7agrcVi0SNHjtRt2rTRe/fu1W3bttUJCQlaa61HjRqlQ0NDdXx8vNZaa5vNprt3764nT56s8/Ly9MCBA3VgYKB+5JFHiuy/sh+KrrV7vk/i/51Ly9YfLN+qp73yoL7wSqTW/wrWqdNv07ajG7S22YwOzyNxtQ9GN4onPhi9JNnZ2WzZsoV+/fpdVqe15siRI7Rs2RKr1Vowv79//36Sk5OJjY0taJuenk5QUBAAM2bMwNfXl7vuuqvIJ4qq4Cnvk6ex2TSmjATOrp5I0K+fE6AspEUNJKT/eLnbpcEq48HowkD+/v4lJnywX3zVsmVLgCK3X2jfvv1lbS8lfKBgSacQFbX7eCpf/7iRYVlLuOHC90TarGS1HQaxzxBSV/64uwJJ+kKIMtlsmnWHEvnuxx/pmfg5r3ttRSsf6Hwv3PQ4tWpHGR2iqABJ+kKIMi1atpS6e6fyntducs21sMaMw9xjHARFlt9Z1DgunfS11gX3lxE1T009XyTKlmnJ54sdJ4j1PUjzP6bzj2MbsPiFYL3xecxdx0CtOuVvRNRYLpv0zWYz2dnZVX4iUly57OxsebKWCzmfYWHe5qMc37qUB2xLaW76CwIjYcCb+HYeBb6BRocoKoHLJv3w8HBOnTpFeHg4QUFBeHt7y1F/DaG1Jjs7m9OnT1OvXj2jwxFOmLjqVxI2L2K0WsE1ppNYghtDnw+gw0jw8St/A8JluGzSDwkJwdfXl6SkJJKTk8nPzzc6JFGIj48P9erVIzg42OhQRCn+SLhI6zAzpv2LGL13AiFep8mt0xr6zML3Wrl61l259Lvq5+dH48aNjQ5DCJehtWbzX8l8uv5XmsYv5tmgNfjlJBLSoBP0moC51a1Qwi22hftw6aQvhHCO1ab5/tcEFqzbxw2JS5jos4pQnwzy6/WE3rOgWW+5VYKHkKQvhBu7tMJNX0wgbeXzfJK/ilo+OVhb3Qo9n8a78Q1GhyiqmSR9IdxQWlYen22LZ+eePXwSvRXvfQsYactDXzccej6JV71rjQ5RGESSvhBu5MyFbOZsOsaOHVu4Xy9jjtcWTHu9oONI1E1PoMIuf9Sm8CyS9IVwE38lZvDMR3N5xGsFL5l2YvP2xxTzCNw4DoIblL8B4RGcTvpKqaHAnY5vv9Jaf1PRnSmlOgI7gWitdXxF+wshitoZn8KxxAzuqnuCFhvfY7nPT9h8Q6DreExdH4aAMKNDFDWMU0lfKTUWeBOYDYQCK5RSo7TW853dkVLKBMx0dp9CiKISL+YwbtEeJv/9evadusCMDUcIPrWOJ32/AX0QFVAXbn4NU8yD4CfXR4iSlZuAlVJBwNvAcK31GkeZBkYDTid94FHgmisJUgiPNiEaMhOpCywG+AAGAH0x4W22YQtuDDdNhOvvAR9/Y2MVNZ4zR92+wLhLCd8hAejo7E6UUg2xf1J4HphSkQCF8HiZiSUWe2ODodMxtRsBXnKPI+Gcci+901qfLzyNo5QKB+4BVlRgP5OB74C4CkcohIfKybMyf2t82Y063i0JX1RIhebXlVLTgL8BPwPvOtlnCNAbaAOUeUtMpdQYYAxAkyZNKhKaEG4jJ8/Koh0nWLp+O0OyV8hZMFGpKvrrtBloCcQC3YCNZTVWSgVin855RmudqJSKKqu91nom9pO9xMTEyM3YhUda/O13hPwynWVeW+0H8fI/QVSiCiV9rfUCYIFS6nPsUzYdy+nyBvCn1nrulYUnhPvLys1nwdbj3GQ6QNtjn3Lf0XVYfWvhFfMQdHsEPmxndIjCjTizescM1NNanyxUHAeMcGL7Q4GmjtU+hR1TSs3TWo9yNlAh3E2mJZ8FW/7ixM8LuNu6gram4xBYD/r9C6+YB8C/tr1hQN2ST+YG1K3egIVbcOZIPxb4UinVXGud4iiLBo470fc2wFzo+wbY/2AMAn6tSKBCuJMvNv3O6bUf83dbHA1VMtm1W0LvKdD+LvD2Ldr42cPGBCnckjNJ/yfgFBCnlHoNiASewb78EqVUa8BS0hW2WuvfC3+vlLrg+OfvWusTVx62EK7nYk4e/jmJ+OycybDts/HVGaTX7wqx0/CPHiD3sRfVotykr7XOU0oNwj6HvxhIAl7QWk91NJkBxAOjqihGIVxaWnYeK39YS+CejxmiNgE2zG2GwI2PE9Sos9HhCQ/j1IlcrfVxYEgpdX2c3Znj04A8qUF4hLTMXFZ//zWRv87gXvZgUb5caHsPYf3+B1WnmdHhCQ8lK4CFqGzWfPhjBUnL3+Ku/MNcNIWS2Olp6vYdh2+tOkZHJzycJH0hKklKaip7V06mT8oSTGnHaRQcRUKHt6jf6wGC5Z44ooaQpC/EVUo5d5I/VrzHtae/oq/KIKXO9dT523/wa30b9eXkrKhhJOkLcYWsiX+y/6s3aJv4Hd3J50DQTaT3f4bGHWKNDk2IUknSF6IitCb7yGb8d07F69B3XIcPO2vfSsPbnqFDqw5GRydEuSTpC+EMm5ULe5aT/tN7NM78Datfbbx6jcd0w2huDJIrY4XrkKQvRFnysknbOo/8zZMJs5wiTddlWYP/ofvwx4gMD8fL6PiEqCBJ+kKUJDMZds7Ctn0mIdnJ7LO1YFXUv+k15EGGhQcZHZ0QV0ySvhCFJR8hY8Mk/H79Am9bDqZWA/k54m6aderPP8ICjI5OiKsmSV8IgFO7yFr3Pn5HvsOsvVime9Dr/teo16IjvYyOTYhKJElfeC6bDf5cRd6mj/A5tY08HcB8PYT09g9yT/+u1AuRC6qE+5GkLzxPXg7s/xLblsmYkg/jFdyI900PkNNuJA/2bU9kiJ/REQpRZSTpC8+RnQo7PyF/63S8s5M46tWC5nfOxnTtUB7VJny9ZS2OcH+S9IX7Sz0O26Zh2z0fU34Wm23t+VQ/RPQNg3iqTWv8vLzwLX8rQrgFSfrCfZ3ZA1smw2/LsaFYlt+dz9TtdOvem4k9mxEWKKleeB5J+sK9aA1//QibP4L4jeR7B+DdfSy2Gx7i3D4Lc25oQp0Ac/nbEcJNSdIXrmVCdCkPCY+A/v+2H9kn/k6qVzjT8kayLWAwK/vfhrdSjO1T7dEKUeNI0heupaSED5CZBMsf4ZRPM97PfZh15l7c2yeaz29qhlLysDYhLpGkL9zGr7FzGLkugP/u24L1N0UR4u9jdEhC1DiS9IXr0LrM6mt73cmWblYCfeXXWojSyP8OUfPl58Lvy2HbtDKbKaUk4QtRDvkfImquzPOwey7smA0ZZyEs2uiIhHB5kvRFzXPuN9j2MXr/YpTVwl7fzgQPeJvm3e4gf0IrvLOTLu8TIA8yEcIZkvRFzWCzweHV9imcYz+TZ/Jlha0X0y39yQtoxWth19LcZML7f//ioc92ERHkx8guTVi44wRJ6TnMuDfG6BEI4RKULufkmFFiYmL0rl27jA5DVDVLOuxdCNunQ8pRdFADPkqP5bPcPnS9tgUjuzTlxhZhmEyy7FIIZyildmutSz0KcvpIXyk1FLjT8e1XWutvnOxXD5gO9Ae8gOXAP7XWmc7uW7ih1HjYPhPbL/Mx5aZzyKcNrUbMRbW5nc5HLzAyMoi6QXK3SyEqm1NJXyk1FngTmA2EAiuUUqO01vOd6L4EaAg8D9QGXgaSgMevJGDhwrSG45uxbf0Ydeg7rCjirF2ZZ7uViBY3MqFlB4K9fOgZHWF0pEK4rXKTvlIqCHgbGK61XuMo08BooMykr5TqB1wHtNVaJzjKIoGhSNL3HPkW+PVr9LZpqLMHyDeHMjP/dtYE3E7/rh2ZHtOYusFyVC9EdXDmSN8XGHcp4TskAB2d6LsL6H4p4TskA3KppCdIP4d15yfkb5+NryWZlFotCLt9ErQdwXUnMnkkOgIvmasXolqVm/S11ucpdESvlAoH7gHmONE3DUgrVtwf2FaxMIVLSdhH5s9T8D24DG+dx3rr9azwG8cNPYZxb+dmmIE+reUh40IYoUJLNpVS04C/AT8D71Z0Z0qpAUBXoG8p9WOAMQBNmjSp6OaFkWxWbAfjMG2fDsc346X8WJAfy8EmI7m5x4180LquHNULUQNUaMmmUuofwP1AF+B2rfXGCvT1B/YDh7XWt5XXXpZsuoicNNK2zIHtMwixJJAf1Ajv7o9wtPFQfIPCaBgqDxcXojpV2pJNAK31AmCBUupzYDLOzetf8i4QBsRWZJ+iZspPPMyZHz4i4sgSQnQ2O2zXsL3eaG6580Fa1a9Nc6MDFEKUyJnVO2agntb6ZKHiOGCEsztRSt0JjANGaK1PVThKUTNojfXIerx2TMfrz9VEahM/mHqS0u5BYmP781idWkZHKIQohzNH+rHAl0qp5lrrFEdZNHDcmR0opToC84ApWuuvryhKYSirJYvDP84hcO9sGuUdg4AIVO//5VD9O7klOhofL5PRIQohnORM0v8JOAXEKaVeAyKBZ7BfbIVSqjVg0VrHF++olPIBvgRSgUVKqcLzTPu11rlXF76oSkln4jn2/Ue0OvkV15DOIaL4rvnL3PL3R/Ey+9PO6ACFEBXmzJLNPKXUIOxz+IuxX037gtZ6qqPJDCAeGFVC93ZAK8e/Nxera+boJ2oQm01jO7UL750zCPt1GWE2K7/4dye/y8N06jmY1j5eRocohLgKTp3I1VofB4aUUtenjH6/ALJOzwUkpmWwZ9V8Ghz6lHa2Q2AOIj/mnyRecz8xLdoaHZ4QopLIrZU9mNaabb/9RcK6mXQ7/zW3qGTOetXnSOeXadF/DGa/YBoZHaQQolJJ0vdAWbn51Eo7gto+nU67FuCLhfiQGM72mEhkzB1gkikcIdyVJH0PYbNpth5JYu+6pbQ/vYieai94+ZLXdjjqpkeIatje6BCFENVAkr6bS8vKY8m2Q6Rt+4whOSt51HSGdHMYmV2eI+DG0QQGhBsdohCiGknSdwOJF3MYt2gPU0ZeT90gP7TWpFvyCc45i94wjRG/fEqIyiK1zrXk9v4XQe3uBG+z0WELIQwgSd+VTYiGzETqYl9Ly3v24nQC+TMwhpjMjYSiyYoeBL3GUbtxV1CymEoITyZJ35VlJpZYHEwG7S2/QPdHoctoaoXKHUuFEHaS9N2U+dmDYJZ71gshipKbprgrSfhCiBJI0hdCCA8iSd9FZeXmGx2CEMIFSdJ3UR9/vRpbaQ89C6hbrbEIIVyHnMh1Qav3H6f/Hy9iMQfj/9g2CGlodEhCCBchR/ouJiEtm4SlL9LedAzvYdMk4QshKkSSvosJPv0zo/iGi9fdj8+1txsdjhDCxUjSdyE6/RwBceMgog3Bd7xjdDhCCBckc/ouYs/xZGyf3831touY7lsBPv5GhySEcEFypO8CMiz5bPr8DTrn7cbS999QT55kJYS4MpL0XcCML5byUO48UpsMwP/GMUaHI4RwYZL0a7hvd/3JsCMvY/ENo/bfZ8hdMoUQV0Xm9GswrTUBa18kynQO/beVUKuO0SEJIVycHOnXYOrXr4nN/oHsbk/i1aKX0eEIIdyAJP0aasP2nehv/wcadSGg/4tGhyOEcBOS9GugnUfOERT3MJZ8GwyfDV4yCyeEqByS9GuYtKw8flv4Ap1Mf8Hgj6B2U6NDEkK4EaeTvlJqqFJqvuPl9PX/SilfpdQ0pVSKUupPpdStVxaq+9NaM3fhfO7L/5rkVnfhd/1/GR2SEMLNOJX0lVJjgbnAOcACrFBK3efkPiYBdwOPAR8CS5RSrSseqvv7Zuuv/P3k66TVakLYiA+NDkcI4YbKnSxWSgUBbwPDtdZrHGUaGA3ML6dvI+CfwH1a6wWOsk7A04BcZVSY1vT/63V8TBmoe5bL4w6FEFXCmSN9X2DcpYTvkOAoL09vwAosLVT2LdDP6Qg9QL7Vhm3HLPyPrsZ7wL/xatjR6JCEEG6q3KSvtT6vtS44oldKhQP3ACuc2H4D4KjWOrtQ2UmgqVLKq6LBuqt5y78j//sXsLW4Gbo9YnQ4Qgg3VqHVO0qpacAhYD/wrhNd/IALxcqyAS8gtITtj1FK7VJK7UpKSqpIaC5r68GT9Nw3Hot3EKZh0+U2C0KIKlXRJZubgd1ALNDNifYW7NM7heU6vl52b2Ct9UytdYzWOiYiIqKCobmelMxczix+ilam0/j810wIdP8xCyGMVaGkr7VeoLUegH1efrITXRKxT/EUVtvxNasi+3Y3Wmu+mD+V4bYfSGr/EH6t+xsdkhDCA5Sb9JVSZqVU42LFccA1Tmx/P9BEKVU48V8P5ACpTkfphpJOH+Uf5yaSFNSWiCFvGB2OEMJDOHOkHwscUEoVvsVjNHDcib57gVPAkwBKKW/sSz1/0lrrioXqRmxW6v4wjmCzJuz+z8HbbHREQggP4cxNXX7CnrjjlFKvAZHAM8DzAI4LrSxa6/jiHbXWNqXUs8AipVQroD7QCehZOeG7npw8K78tepHOJ7aghs1AhbcwOiQhhAdxZslmHjAISAIWAy8DL2itpzqazABeLaP/YuB27Kt1soABWustVxW1C1u4ZDEdj0wnqdkd0OHvRocjhPAwTt2+UWt9HBhSSl0fJ/p/B3xXocjc0Mb9hxlw8CXSfOsT8bcpRocjhPBAcs/eapJ0MYecZY8RqVKxjlwFfsFGhySE8EBya+Vq8t38d+ivt5LabTy+UV2MDkcI4aEk6VeHpEPcmzqNpIjuRAwYb3Q0QggPJkm/iqVnpMOS/8bkG0DEfZ+CSX7kQgjjSAaqQtm5VtZMegTOHYChH0NQpNEhCSE8nCT9KrR40SzuzP2GM63vh1a3GB2OEEJI0q8q63buZ/DRNzhXK5oGI5y5IakQQlQ9WbJZBRIuZOIfN5ZAlYu67zPw8TM6JCGEAORIv0rkrP+AbhzgYuwbmCPbGB2OEEIUkKRf2U7tptn+D7C2uYOIXqONjkYIIYqQpF+JDhw5SfqC+9BBkXgNmSRPwRJC1DiS9CtJhiWfswsfpVb2GbJvnwH+oUaHJIQQl5GkX0lWzn+f/tYNJHR8glotexgdjhBClEiSfiVYu3krQ069x6ng62k05GWjwxFCiFJJ0r9KFzMzqbfmUbTJm8hR88HkZXRIQghRKkn6Vyl489tcxxEst36Id50mRocjhBBlkqR/FVL3r4ItkyDmQcK7/JfR4QghRLkk6V+hvX/8Sf7XY8gIbgkD3jQ6HCGEcIok/SuQlmUh+6uHCFZZqBFzwFzL6JCEEMIpkvQrSGvNmrmv0d32C4ndXiagSQejQxJCCKdJ0q+gtet+ZEjidI6F9abxLY8bHY4QQlSIJP2KsGRww66nyfAOpckDc+Q2C0IIlyO3Vq6IVf9LSNYJbPeuwBQYbnQ0QghRYXKk76QfFk+DPZ9Dz6cxtehtdDhCCHFFnEr6Sql6SqllSqkMpVS2UmqRUirAyb5RSqk1SqmLSqlzSqlpSimXeqrI7r176Pbbvznufy30ec7ocIQQ4oo5O72zBGgIPA/UBl4GkoAyz2QqpRTwBXACuB1oAkwA8oAnrizk6pWanoXPiocwKUXdBz4DLx+jQxJCiCtWbtJXSvUDrgPaaq0THGWRwFDKSfpAc6ArcLvWOsnRty4wDhdI+lprtsx5lkH6ECf7TqFx3RZGhySEEFfFmemdXUD3SwnfIRlw5pD30tnOwvvxAyzOhWes/Zu+5daUBRyqfweNe91rdDhCCHHVyk36Wus0rfXBYsX9gW1ObH8/kAi8pZQKVUpdD4wFFlY40uqWlUL7HePJDo4i+v6pRkcjhBCVosJLNpVSA7BP2fQtr63WOlspNQJYDzzgKI4DSrxZjVJqDDAGoEkT4+5YmZObj+2rh6mVmUTAP38EvyDDYhFCiMpUoSWbSil/YCrwvdZ6nRPtfR3tNwAjgWeA7sCHJbXXWs/UWsdorWMiIiIqElql+umzt6h1bDXpPV+CBh0Ni0MIISpbRY/03wXCgFgn2w8B6gJdtNY5AEqp34E4pdRbWuszFdx/lduxfSP9TnzE4eBuRPeW2ywIIdyL00f6Sqk7sa+6Ga21PuVkt5bAsUsJ3+EAoICmTkdZTZJSLhD2/SNkmgJo/OCnYJJr14QQ7sXZi7M6AvOAKVrrryuw/SSgTbELuQY6viaU0N4wWmv2zxlHC06SedtU/GrXNzokIYSodM6s0/cBvgRSgUVKqZhC1fuBBoCv1vpQCd1/xL5Ec7NSapWj7V3Aaq11/FXGXqnyf/uGfhnf8FvU/Vx7w2CjwxFCiCrhzJx+O6CV49+bi9U1A14FooA+xTtqreOVUn2B94H/wX4l7mrg4SsJtsqkncLn28egfkeuvWei0dEIIUSVKTfpa61/wT4HX5pR5fTfAnSrWFjVJzsnlxMf303L/Fy8RswBb7PRIQkhRJXx+DOVmz99ntY5+zlyw6sQJrdZEEK4N49O+ts3xNEnYQ6/hQ2g1YAxRocjhBBVzmOT/rlz52i87nHOe0UQ/eAseQqWEMIjeGbS15rkLx4mQqeSP2w25oBQoyMSQohq4ZlJf89ntE39ieQuz9CoXS+joxFCiGrjcUk//uAe9Hf/C816EXmrPAVLCOFZPCrpZ2ZmkL/4AdKt3uhhM+Q2C0IIj+NRWe+XOU/Q0naMM33eQwU3MDocIYSodh6T9HesXkjP5CX8EnkX1/S+y+hwhBDCEB6R9M+cOkbLreM55tWMdg98ZHQ4QghhmAo/Ocvl2GxErHkCrXLJ/vun+PjWMjoiIYQwjNsnfb1lEj7HN8DtH9EwuqPR4QghhKHcenrnj13rsP74b7JbDoJO9xsdjhBCGM5tk/7FtBSC4h4mWdUmf/BHcpsFIYTA3aZ3JkRDZiIAwY4XALNuhGcPGxWVEELUGO51pO9I+E6XCyGEh3GvpC+EEKJMkvSFEMKDSNIXQggPIklfCCE8iHsl/YC6FSsXQggP415LNmVZphBClMm9jvSFEEKUSZK+EEJ4EEn6QgjhQSTpCyGEB5GkL4QQHkRprY2OoURKqSTg+FVsIhw4X0nhuAJPGy/ImD2FjLlimmqtI0qrrLFJ/2oppXZprWOMjqO6eNp4QcbsKWTMlUumd4QQwoNI0hdCCA/izkl/ptEBVDNPGy/ImD2FjLkSue2cvhBCiMu585G+EEKIYiTpCyGEB3GJpK+UekAptb5YWTOl1NdKqXSl1Fml1GtKKVNl1RtFKVVPKbVMKZWhlMpWSi1SSgU46pQjzrNKqRNKqQeK9b2qeqNU5ZgLteunlDpSHeNxRhW/z4FKqblKqQtKqVyl1GqllKH3F6/q91gp1VspNVsptVApdV91jass1fF77WjbQCmVppTq41RgWusa/QK6AlnA+kJlAcAxYDXQB3gQuAhMrIx6g8e7ETgKPAa8AuQBkxx1zwMW4FHgbiAdiC3U96rq3XHMjjYtgXNAvNFjrab3+TMgERgPPOX43V7pxuO9HcgEpgAfAjnAK+78Hhfbz9eABvo4FZfRP5hyfmh9gTTgF4om/X86yoMLlT0FZAPmq603cLz9gFSgfqGyacAJwNfxn/fFQnUvAT84/n1V9e44Zsf37bAn/F3UkKRfxe9zSyAXuK5Q/XhHwvF1w/GagCPA6EL1LwAn3fU9Lraf2x1t3Sbp/wu4A3iVokl/GrC2hMFroMHV1hs43hDgmmJlrwMJwE2O+FoWquuI/ajG+2rr3XHMju8fcrxGUXOSflW+z7WAdsW2/d+ADfB30/E+CPgUqn8ASHTX97hQWQD2W9U8SgWSvuFz2OV4XWu9ooRyKxBWrKwN9qOZpEqoN4TWOk1rfbBYcX9gG/Y/VjnYj2ouOYn9qKBRJdQboorHDDBLaz2jCkK/YlU5Zq11ltb6QAnb3qe1zq68UTivGsY7R2udB+CYMx8NlJQ3qk01/F6D/Y/IaeDjisRWo5O+1tpWStVGoINS6r8BlFJtgGeA1Y43/2rrawSl1ADs5zQmAX5Amnb8iXe49J84rBLqa4RKHnNZv0M1RmWPudi22wAjHNuuEapqvEqpV4DDgMI+XVtjVPaYlVKdgIeBf1b0d7xGJ/0yLHW8ZiulUoHfgAhgeiXVG04p5Q9MBb7XWq/DflLHWqxZruOrfyXUG64KxlzjVeWYlVIK+5Wdv2E/uWu4Kn6PdwObsU+F3FpJIV+1yh6zUsoL+/s6QWv9e0XjcckHo2ut84HhSqkbgVbAm8BZrXVcZdTXEO9i/6se6/g+EairlPLSWl/6hajt+JpVCfU1QWWP2RVU5Zifwn502cXxO18TVNl4Hf9/45RSbwDTlFJLa8i4K3vMj2E/l/HmFUVj5MkOZ18UO5FbrC4G+0mMwVVRb9B473TENLxQWTj2k3HdCpUNcLSLvNp6dxxzse2PooacyK2OMWNP9rnA00aPs4p/r72AqGL76e7Ov9fAese/S3qtLzcmo38oTv7gXi1tMMBiYHMZfa+q3oCxdsS+JndyCXUbga8Kfb8M+K2y6t1xzIXKR1GDkn4Vv8+NgTPAShz31zL6VVXjBVoD+UB0ofr7sC9jNGz5dRWPuaVj24VfGvtS9JblxmX0L4OTP7xXKSHpY1+DbQV6ltLvquoNGKcPcAj7Wt4bsX8KufQyY1/KlQv8CKx1vNF/L9T/qurdccyF2o2ihiT9anif12G/DuXmYtsOctPxrgYOAkOB4dj/4E1w5/e4hP1p3GGdfqEBvUrJSX8p8G0Z/a6q3oBxdqL0j21RjjZdgVXAFuBvJWzjqurdccyONqOoOUm/ysYM1Clj233cbbyOujDsJ6pTgVOOfGHYtSfV+XtdqK3T76/cWlkIITyIqy7ZFEIIcQUk6QshhAeRpC+EEB5Ekr4QQngQSfpCCOFBJOkLIYQHkaQvhBAe5P8AHekfPM4FG4AAAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#15.1\n",
    "\n",
    "import numpy as np        #GM(1,1)前面的1是表示一阶微分方程，后面的1表示一个未知数\n",
    "import sympy as sp        #该GM（1，1）模型只能描述单调变化过程\n",
    "from matplotlib.pyplot import plot,show,rc,legend,xticks     #导入数据库\n",
    "rc('font',size=16)\n",
    "rc('font',family='SimHei')      #使plot画的图中可以显示中文\n",
    "x0=np.array([2.67,3.12,3.25,3.41,3.56,3.72])     #导入数据\n",
    "n=len(x0)\n",
    "jibi=x0[:-1]/x0[1:]      #求级比\n",
    "bd1=[jibi.min(),jibi.max()]       #求级比范围\n",
    "bd2=[np.exp(-2/(n+1)),np.exp(2/(n+1))]      #q求级比范围的容许范围\n",
    "x1=np.cumsum(x0)      #求累加序列\n",
    "z=(x1[:-1]+x1[1:])/2.0\n",
    "B=np.vstack([-z,np.ones(n-1)]).T\n",
    "u=np.linalg.pinv(B)@x0[1:]     #最小二乘法拟合参数\n",
    "sp.var('t')\n",
    "sp.var('x',cls=sp.Function)      #定义符号变量和函数\n",
    "eq=x(t).diff(t)+u[0]*x(t)-u[1]     #定义符号微分方程\n",
    "xt=sp.dsolve(eq,ics={x(0):x0[0]})      #求解符号微分方程\n",
    "xt=xt.args[1]     #提取方程中的符号解\n",
    "xt=sp.lambdify(t,xt,'numpy')      #转换成匿名函数\n",
    "t=np.arange(n+1)\n",
    "xt1=xt(t)     #求模型的预测值\n",
    "x0_pred=np.hstack([x0[0],np.diff(xt1)])     #还原数据\n",
    "x2002=x0_pred[-1]      #提取2002年的预测值\n",
    "cha=x0-x0_pred[:1]\n",
    "delta=np.abs(cha/x0)*100\n",
    "print('1999~2005的预测值：',x0_pred)\n",
    "print('\\n-------------------------------\\n','相对误差',delta)\n",
    "t0=np.arange(1999,2005)\n",
    "plot(t0,x0,\"*--\")      #画出原来的数据图\n",
    "plot(t0,x0_pred[:-1],'s-')      #画出由灰色预测模型下预测的数据图\n",
    "legend(('实际值','预测值'))    #显示两幅图的图标\n",
    "xticks(np.arange(1999,2005))\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "********** 级比检验\n",
      "[0.73676681 1.46848739 0.7020649  0.98083183 1.78387097 0.55258467\n",
      " 1.87       0.47468354 1.17037037 1.32939439 1.29445507 0.54479167\n",
      " 0.98025868 1.84489796 0.79565326 0.6777853 ]\n",
      "[0.47468354430379744, 1.87] [0.8948393168143698, 1.1175190687418637]\n",
      "********** 建模\n",
      "[ 691.6  1161.6  1623.2  2170.9  2602.4  3037.9  3468.4  3934.4  4520.4\n",
      " 4993.5  5353.5  5798.4  6380.2  6833.25 7192.65 7688.1 ]\n",
      "[[-6.91600e+02  1.00000e+00]\n",
      " [-1.16160e+03  1.00000e+00]\n",
      " [-1.62320e+03  1.00000e+00]\n",
      " [-2.17090e+03  1.00000e+00]\n",
      " [-2.60240e+03  1.00000e+00]\n",
      " [-3.03790e+03  1.00000e+00]\n",
      " [-3.46840e+03  1.00000e+00]\n",
      " [-3.93440e+03  1.00000e+00]\n",
      " [-4.52040e+03  1.00000e+00]\n",
      " [-4.99350e+03  1.00000e+00]\n",
      " [-5.35350e+03  1.00000e+00]\n",
      " [-5.79840e+03  1.00000e+00]\n",
      " [-6.38020e+03  1.00000e+00]\n",
      " [-6.83325e+03  1.00000e+00]\n",
      " [-7.19265e+03  1.00000e+00]\n",
      " [-7.68810e+03  1.00000e+00]]\n",
      "[1.83087632e-03 4.80930834e+02]\n",
      "********** 解方程\n",
      "Eq(x(t), 262677.947893454 - 262265.947893454*exp(-0.00183087631680452*t))\n",
      "********** 求预测值\n",
      "1989~2006的预测值： [412.         479.73720894 478.85967303 477.9837423  477.10941382\n",
      " 476.23668467 475.36555191 474.49601263 473.62806392 472.76170285\n",
      " 471.89692654 471.03373207 470.17211657 469.31207713 468.45361087\n",
      " 467.59671493 466.74138641 465.88762247]\n",
      "\n",
      "-------------------\n",
      " 相对误差 [ 0.         14.21008424 25.75096455 11.87615371 13.72343331 53.62473699\n",
      " 15.2646075  58.16533754 25.05885065 12.45153651 16.17354174 50.10635184\n",
      " 18.37289643 20.13068803 47.08119651 16.81157005 20.97165824]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABoyklEQVR4nO2dZ3gc1dWA39ld7apXS1azLffe5QIY4wI2GDCmxIEEEkPAdFIIhBZCwIGPhBBiWmxaQgBDKKYZm2JsXHHvvav3Lq22zvdjdqXVanuXNO/z7CNp7p0zd1Yz57ZTBFEUkZGRkZHpGSjC3QAZGRkZmdAhK30ZGRmZHoSs9GVkZGR6ELLSl5GRkelByEpfRkZGpgehCncDnNGrVy8xLy8v3M2QkZGR6VLs2rWrShTFdGflEav08/Ly2LlzZ7ibISMjI9OlEAThnKtyeXlHRkZGpgchK30ZGRmZHoSs9GVkZGR6ELLSl5GRkelByEpfRkZGpgchK30ZGS+paGhl4bKtVDS2hrspMjJeIyt9GRkvWbr2BDvO1rD0uxPhboqMjNdErJ2+jEykMfSx1eiM5ra/39lWwDvbCtCoFBxbclkYWyYj4znySF9GxkM2PjiTeaMy2/6OjlJw1bhsNv5hZhhbJSPjHbLSl5HxkIzEaJQKAQCVQkBnNJOgUZGREB3mlsnIeI68vCMj4wVlDdLm7aT+KSRGR1HZpAtzi2RkvEMe6cvIeMFN5+UBsPVUDRqVkmU35Ye3QTIyXiIrfRkZL6iwjPQvH5PFt4fLadEbw9wiGRnvkJW+jIwXZCRGM31IOjdO6YfWYOK7IxXhbpKMjFfIa/oyMl4wf2w288dmYzaL9E7U8PneEuaPzQ53s2RkPEYe6cvI+IBCITB/bDZagxGzWQx3c2RkPEYe6cvIeMGClzczKieRJQtG8/Blw1FYTDhlZLoK8khfRsYLzlU3t/1uVfjNOnkzV6brICt9GRkP0RvN1LYYOjhjfbK7iAlPfSsHX5PpMshKX0bGQ6osjljpCZq2Y6NzktAZzXy1vzRczZKR8QpZ6cvIeEhFo6T0M2yU/uDeCQzLTODzfSXhapaMjFfISl9GxkPiNUqunZBL/15xHY7PH5fN7oI6CmtawtQyGRnP8UrpC4IQLQjCSUEQnrE5JgiC8GdBEMoEQSgQBOFmu3NclsvIdBUGZSTw94VjGZAe3+H4lWMkO315tC/TFfDWZPMPQBzwtM2xhyyf3wE1wHJBEM6KorjOw3IZmS6B3mgmSikgCB3NNPukxvLstaM5f2CvMLVMRsZzPFb6giD0R1Le94ii2Gg5pgEeBp4URfFly7GBlmPr3JUH8kZkZILNnz4/xIbjlWx+aFansp9O6huGFsnIeI83yzv/BI4Cb9kcywcSgA9sjn0JTBcEQeVBuYxMl6GyUUdCtPPHds3BMr6Ql3hkIhyPlL4gCJcCVwJVwDuCIPxREIQ4IBtoBU7ZVC8ENECuB+UyMl2GysZWMhKdJ0z5749nee6bY4iiHJZBJnLxdKRv3bjtDWQBfwI2A9FAvdjxKddafqZ5UN4BQRAWC4KwUxCEnZWVlR42TUYmNFQ06kiP1zgtnz82m3PVLewvqg9hq2RkvMOt0hcEYRIwDnhNFMUxoijOBGYBo4D+gMnuFL3lZwygc1PeAVEUl4uimC+KYn56errHNyEjE2zMZpGqJh0Zic6V/qUjs4hSCrIVj0xE48lIf7Dl5/PWA6IobgBOADcAGYIgKG3qp1h+tgAVbsplZLoERrPI7dMHcoELC52k2CguGpLBl/tL5MibMhGLJ0q/yfLzjN3xVuB7QAlMsjk+3vKzBDjoplxGpkugVin4/dyhTBvs2ixz/rhsYqKUlDbIsXhkfKOioZWFy7YGLZ6TJ0p/JyACY60HBEHIAIYCW5HW9u+3qX8ncFgUxTJRFKtclfvZdhmZkNGiN1LbrHe7SXv56CzW/X4GOcmdVi9lZDxi6doT7Dhbw9LvTgRFvluzSVEUSwRBWAH8VxCEh5DW6f8IlAMfIc0A1gmC8B0gIK3332Aj4iE35TIyEc/qA2Xc/+E+1v9+Bnl2YRhsUVrCLeuNZhQCqJRypJNAUtHQyj0r9vDSz8Z3iHbaHRj62Gp0RnPb3+9sK+CdbQVoVAqOLbksYNfx9Im8BfgMWAp8iLQZe5koiq2iKG4GLgSMSJuz14ui+L71RHflMjJdgUoHETadcaS0gfwl37LhhGyBFmiCPQoOJxsfnNkh9WZ0lIKrxmWz8Q8zA3odjxykRFHUAQ9aPo7KtwGXujjfZbmMTKRT0aAjTq0kTuP+lRmYHo8gCHyxr5RZw3qHoHXdG5NZZNgfV2MwtS+tBWsUHE4yEqNRKaWZokohoDOaSdCoAj6jkb1iZWQ8oMKNY5YtapWCy0Zl8sW+ErR6EzFqpfuTZDpxuKSBlXuK+GxvCQaTyMXDM1h3tAKTKI2C547M5NHLh4e7mQGlWWfkxil9uX5yX97fUUhlEDZzZaUvI+MBlW4cs+y5cmw27+8o5PujFVw+JiuILet+HCyu5/cf7uNoWSNRSoEZQzO4enwOG45XYhKljcFgjYLDzbKb8tt+X5KTFJRryEpfRsYDbpzaD8GLHOhTB6SRnqDh833FstJ3Q7POyJqDZaTFq5kxNIOspGjiNSqeumokl4/JJjVODcBne4sZ0juekrpWFozPCcooONw8/80xGlqNPDF/ZNCuISt9GRkPuNJmg80TlAqBP88fSW8Pl4R6GkaTmU0nq1i5p5hvDpWjNZi4cmw2M4ZmkBav4aM7z+90zrKb8nl1/SmeXXOUhy4bRrwH+ytdje+PVZASqw7qNbrftyYjE2AMJjMnypvomxbrlaKZN1oe4YNjM8tb397J+mOVJMVEcfWEHK4Zn8PEfiluJEF2snR+aZ2Wwb0TgtruUGM2i5ysaOJnk/sF9TqyEbGMjBtK6rTMW7qRNQe99yfcX1THf7acDXyjuhDPrjnK9jM1XPrCRhpaDQD88rw8/nXjBLY/Opunrx5Nfl5qp+Q0jhiUEc/sYRnBbnJYKKxtodVgZkjvePeV/UAe6cvIuMGaEN0TG3171hwsY9mG01wxJos0LzaCuwP2zkY1zXrGPPGNX2aWI7OTeGPRJPcVuyDHy6WIN8GewcgjfRkZN1RalH6GD0p//rhsTGaRrw6UBrpZEc/GB2cyb1Rm29+BdDbqjjkLRFFkWGYCg4M80peVvoyMGyoswdN8UfrDMhMZ0ju+R4ZbzkiMbrN4CqSz0TWvbOb+D/cFoIWRxZyRmaz5zXQSo6OCeh1Z6cvIuKGySYdKIfhsVTF/bDY7ztZSXKd1X7mbUddioE9KDEtvGMfPp/RrC2fhDyqlgqKanvddBgp5TV9Gxg3zRmcxKCMehcILQ30brhybzWsbz3CivLHHRd9897apbb/PG+2d2aszcpJj2H6mJiCyIgWTWWTas99z54yB/OK8vKBeS1b6MjJuGJmdxMhs370j+6XFseuxi3tkxE2zWfS5s3RGTnIMZQ2tGE3mbvOdFtS0UFrfSrQq+CE7usc3JiMTRLadruZsVbNfMlRKBaIoojPaZw/t3tz7/h6ueWVzQGVmJ8dgMottVlXdgePljQBB38QFWenLRDDBziDkKXe/t4dlG075JaPVYGL233/glXX+yelqFNVqiVUHdkFhdE4SN07ti8KbuBgRzsmK0Jhrgqz03RIpiqcn8szqI2GPnW40malu1pHup8VJdJSSjEQNX+wv6Zbmhs4ortWSmxLYfYzRuUksWTCazKTuE+LiuGW/JxShJWSl74a/fXMs7IqnpzH0sdXkPbSKlXtKEEUpdnreQ6sY+tjqkLelplmPKPrmmGXP/LE5nK5s5lBJQwBaFvm0GkxUNemCsnltNJlp1hkDLjdcjM1N5tqJuSG5lqz0nWBVPB/uLAq74ulphCqDkCdU+OGYZc9lozJRKQS+6CE2+1YT1ZwAj/QBpj6zlqe/OhJwueHilmn9+d0lQ0JyLVnpO2HjgzMZmZ3Y9nc4FU+wibQlrIzEaOI0khWDQghv7HTrdxKIkX5KnJoLB/fii30lmM3df4knOkrJovPzGBWEuPC9E6Mp6SZ+D60GE62G0G3wy0rfCRmJ0cRZNqDCrXiCTSTmHd11rhaQoioGyqnHF8b1SeHfN09icEZgrCruuGggf7xiBN1f5UumlU/MH8mQIGxO5iTHdBtnt/XHKhn++BoOh2jZT7bTd0FKXBT90+LISYkmr1d8t0vaYB8QK5Lyjk7un8rx8iaadSaWLBgVtnakxkmJPQLFlAFpAZMV6dS3GIhWK9AEwfY8OzmGzSerEEXRo+ickcyJ8kZEEfqlxYbkerLSd4LZLPKvGyd2+QfKEWermll1oJR+abEcL29CqRAwmcWIyjtaaHGzT4iOotVgIjoqPHlmd5ytoVlnDKjiL6nT8tneEm67sH+3cS5yxBNfHGLH2Ro2/WFWwGXnpsTQrDfR0GokKSa4sWqCzfGKJnJTYogLUVKY7vvE+cnmU1VMfnptyKZcoeKmN7Yx47n1/O3rY8RpVEzom4zZLKJRKSJqCauotoXLRmWy4cGZYVP4AG9sPMOSVYHdMDxQXM+za46y+VR1QOVGGsW1WrKDFHZi6oA0Hpg7NCiyQ82J8sagLIE5Qx7pO2FPQR1VTTriNSqueWUzvzgvjwXjc8LdLK84VdnEV/tL2VNYxxu/zEcQBC4c3IuLhqRz2egscpJjuP2/O7lwcC8m9U+lvEEXEUtYoihSVKtlVgQky6hobA2I5Y4tM4amkxCt4vO9JVw0JD2gsiOJ4jotU/qnBkX2qJykoGwQhxqjyczpymYuGhq650BW+k7YW1jHoPR4clNiOFLayL6iui6h9ItqW/hkdzFfHSjlaJnk2j2xXwp1LQZS4tQsnj6wQ/1lN+Vzy7938NWBMlb/+sJwNLkTepOZX03rz+icJH7x5nZ+MjHX6xy1gaKyScfEvu7T+HmDRqVk7shMvj5YRqthVFhnMsHCYDJTWq8NirkmSAOD0vpWVEohImamvmI0izx02TDG9kkO2TXl5R0HiKLInoJaxvVJRqEQGJqZwNHSxnA3qwO2ZpYnyhvbYr4fKmngH98dJyFaxZ+uHMGPD8/m4zvPJyXOeVjgvqmxFNa0RIynqEal5MFLhzF3ZCY/nqoOmzOTKIpUNOgCYq5pz/yx2TTqjKw/VhFw2ZFAWX0rZpGAe+NaEUWY8dx63th4JijyQ0V0lJJbpvX3KD9woJBH+g44V91CbYuB8ZYR3vCsBNYcLIsoS4Elqw6z/UwNF//9Bxpajfzm4sH85uIhzBiazo8Pz6Z3ouejn9yUGJp0xrbZQLix5lFNjI4iKzmaotqWMLXDiM5oDspI8vyBaWQlRXOgqIE3N5/tkDS8OxCnUfHovOFM7Bec5R2FQiA7KbrLm22eqWpGpRDokxoayx2QR/oOUSoEbr4gj6kDpAd2aO8EalsMERHVz+op/Pk+Kf1eQ6vkim4N5KVRKb1S+CCN9EEK7xoJ/HfrOcY88Q1avYnclPDZY8eplXz9m+lcNT7wS0sqpYIfHphJvVYfcT4SgSA1Ts1t0wcwKED+DY7ICeOzESie++YYP399W0ivKY/0HdAnNZY/XTmy7e9xfVO4ZERvtPrwh8Xd+OBM7n5vNzvOSs5LgTCz7GuxDy6q1YZ0bdEZhTUt9IpXE6NWkpMcw7pjlWFph0qpYGhmcKwqItlHIhAU1bYgigR1BJudFMOGE+F5NgKFZLkT/HDKtsgjfQecKG/EYGp/Icf1Sea1X+ST1ysujK2SyEiMZkjvBAQImJnloPR49j5+CfNGZ7qvHAIKa1vITZGUxcjsJAamx2EKQ9iCI6UNvL31LE1BCOy18cGZzB+XjTW/SHcL87F07QmueXVLUK+RkxJDRaMOvU3n2ZUwmMycqWoOSThlW2Slb0erwcS8pRt5/tvjncoiJQFGVZOOn0/tx8q7LghIiAKVUkFyrDpi9iuKbMLx/vL8PN5ffB7KAGdf8oTNJ6t4/LNDmEyB73AyEqNJ0Kiw9mWR5CMRCIrrtEFPDTlnRCbPLxyL2EWDWpytasZgEuWRfrg5VNKAwSQyzm6Z4/7/7eOqlwKbAcgXRFEkNU7DZaMyGZGdyJIFo1h2U77fcldsL+DldScD0EL/MJlFSuq0Id3YckZFow61SkFiTHBWQauadORbrDauHpcTtvhCwaC4NnjmmlZGZCdy9fjcoIR5CAXHyy2JUzLkkX5Y2VMgrZWPt1P6GYkaTlU2hX0qWVSrZcX2As74mb7Pnq2nqvlgR2FAZfqC0Wzm8StHcsmI3gBUN+m45Pkf+HhXUcjbUtEgOWYFawa07KZ8npgv7R1dNDQ9IJ13JGA2i5TUtQbNXNOK0WRm17kaCiPEAMFbJvdP5ZWfTwjqZrcjZKVvx97COnKSY8iws4AZlpmAwSRyuqopTC2TOFhcD0gp4wJJ39RYiuu0GE3h7dQ0KiU3Te3HBIu5bFJMFKermgPeyXlCZZMu4N649gzPSiReo2L7mZqgXieUVDbp0JvM5AZ5ecdoFrn21a2s3FMc1OsEi/QEDfNGZ4XcOU9W+nbsKajrtLQD0ssJhN1J60BxPSqLw1gg6ZMqJZsurQ9vGIbSei3HyhrbNm5VSgWZieGxxw6WY5YtSoXA7dMHMDlI4QrCQbxGxUs/G8+Fg4MbWiA6Skl6gqbLxtX/bG8xx8pCr09kk00bRFHkr9eNIVbdueft3yuOKKXAkbIGFhC+cAwHSxoY3Dsh4KMD6xp6YU1LWNfTP9hRyD/XnuDoU5eiVEj3mJMSQ3Ft6F/sj+48PyTLeffOHhz0a4SSOI2KK8aEJmxGdheNq683mrn/f/u4/aIBPJA5LKTXlpW+DYIgcMGgXg7LopQK7ps1OOxBnnQGE+P6BL4NfVNjiVUrqdMaAi7bGwprtPROiO6wOZebHMO2MCx/hDJkb0VDKyJ47VgXiRwra6S2Rc+U/qlBtwjLTY7hSFnXi4R7pqoZo1kMaXRNK7LSt2HTiSrMosh0J5EPI2FE9sHt5wUlRk5OcgyH/jw37GabRbUt9EntuBY8dUAaCoUQ0jAY1U06Xt90hqvH5wT9xdQZTUz76zoWnZ/HI/PCn8vAX97eepZVB0rZ+/icoF8rOzmatUfLIypEiiccL5eWdUJtuQOy0u/Ay+tO0qw3OlX6RpOZs9Ut5CTHEONgCShUBOPhjpQXpqhW22l9e+GkPiyc1Cek7SioaeHV9afI75cSdKWvUSkZm5vUbTZzQ2Gjb+Wnk/oye3hvRBEi5BH2iBPljSgEGJAeeodPeSPXgskssq/I8Saulc2nqrn4+R/YW1gXsnbZsnzDKRa9tT1oSbWXbzjFIysPBEW2J1jD8fZxYOonimJIvXKtcZZC5Sw1KS+Vg8X1tOgD7/0baoprQ6f0B2XEt80EuxLHy5volxYXlrDastK3cLy8kRa9ifF9k53WGW6xmDkapjXEraeqKatvDdoDfqaqhTUHy4Ii21OW35TfKXZ+cZ2W0U98E1LTvDalnxhc6x0rk/qnYjSL7C2oC8n1goU1AU6wHbOstOiNfLGvhNOV4TWl9pa//mQMr/8yPH4ZstK3YB29j+/jPK51eoKG1Dh1WMw2RVHkQHEDI7ODt5HcNzWWmmZ9UGLNeEKUUsHFI3p3ikXSK15Nk84YUgueykYdggBpIQo1PbFfCoIA28927SWe2hYDWoMpZCN9ncHMvSv2hC0on68kRkcxMD20TllWZKVvYX9RPSmxUS4z0guCwLDMBI6Wh17plzfoqGrSMTonMWjXsG6ghsvD8WRFI+uOVXQIdgfSmndGgiakcfVrm/WkxalDlrg8MTqKf904kYX5od27CDTxGhUf33k+80ZnheR6ybFRxKqVYTHp9ZXCmhae+/pY2PJEyErfwpIFo/j8nmluNzSHZSZy3MZ5KFRYPXGDaTIa7rj6n+4p4db/7MTRfyDUsdOfWjCKjQ/OCtn1AOaOzAxaIvFQoVYpmNgvJWT3IQgC2ckxXcpBa3dBLS+tOxm2GbWs9C0oPcxec93EXJbeMD7kqQU1UQouGJTGiOzgjfT7psYyMD0ubGkTi2pbyEqKdji6zgmDE06oLbTqtQbe3Xauy61P27KnoJaVe4pCOigKx7PhDyfKm1AqBPqHKVS7rPSBfYV1PPbpgbY8s64YkZ3IJSN6h2zab+XCwem8e+tUYtXBs7JNjlWz9v4ZXDoqNFNzewptQirbM2dkJgvGhc4T+onPD/Hl/pKQXQ8kL81HVx7kuyPlIb1uIPl8XwmPrTxIKI1putpI/3h5I/3SYsMWHVRW+sCmk1W882MBGg/Np7acqmLXudogt6ojrYbIiOUfTIpqW+iT4ni2NX9sNr+9ZEhI2mE2i7zz47mQJ2RPT9DQv1cc28+E9tkKJNaQyqH0+7hn1iBW3nVByK7nLycqmhgSBqcsK7LSR5qSDkyP89jt/o+fHmTZD6eC3Kp2KhpbGfmnr/koBOGFn//2ODe9EdqcnSB1auUNuraMWY6o1xpCkrKytkWP0SwGPcKmIyblpbDzXE3QfDGCTVEIbfSt5CTHtKX8jHQMJjPVTbqQJ06xxa3SFwThekEQRLvPd5YyQRCEPwuCUCYIQoEgCDfbneuyPBIQRZE9BXWM7+vcVNOeYZmJHA1hdLxDxQ2YzKJDp6VAo9Ub2X4m9EonSqngq/su5Cf5uQ7LT1Y0MvbP3/DN4eD7EYTaMcuWSXmp1LUYOFHRNdf1i+tCZ6NvpaZZz2sbTnOyC3xnUUoF+/40h7tmDgpbGzxZIB4HfA48ZXPMOu99yPL5HVADLBcE4awoius8LA87hTVaqpv1Lj1x7RmWmcCqA6U06YzEa4IfyeKAxXJnZAiCvfVNjUVnNFPZpAtp8C+lQnC5SW21BikKgWlepUXpBzussiOm9E8DJAfAYCVlDxZNOiP1WoPL2VowaNEb+ctXR0iMUYU8IYkvCIIQFk9cK54q/e9FUdxpe1AQBA3wMPCkKIovW44NtBxb5648YHfgJ+WNreQkx7j0xLVnmCW2/rGyRib283yG4CsHiusZ0CsuJB1Mrk2I5VAq/T0FtRwpbeS6ibmoVZ0noLFqFalx6pBYabToTZZ8taFX+n1SY9jzx0tICZFTWCCJjVKy5aFZIVdovROjUQh0CVt9aa+onqevHh22eFeerOmPB3Y5OJ4PJAAf2Bz7EpguCILKg/KIYFJeKpsfmsWILM9NIYdZRmChSoBwqLg+ZCGdw2Wrv+ZQGX/6/CAqF2YfOcmhiat/6ahMDvx5LnlhMKkTBKFLKnwAhUKymU8NcfujlAp6J0ZTXBfeBECesP5YJTvO1oY1wKFLpS8IQjaQAfxFEIQmQRBKBEH4i0VpZwOtgO2OZiGgAXI9KI8ovPkn5KbE8Pk9F3DNhOCbEJrNIrdM68/V40NjrpiTHMO0Qb1CGkse2jcAXcUVyk2JCZsXYyg5WFzPzW9t73L3+uPpal5edxKdMfSWZlIylcj/vk5UNIZ1Exfcj/QnAmak5ZirgOeA3wN/AKKBerGjJ491GJbmQXknBEFYLAjCTkEQdlZWBj+Whs5oYsbf1vHJbu+sYgRBYExuckimsQqFwK0XDmDmsIygXwukFHTv3DqF2cN7h+R6Voo8yNj1k/xcFk8fEPS2LF17gr+sOhz06zhDIQisO1bZ5UItrztawT+/O0GUIvRGgTnJMZRE+EhfqzdRUNMSlhj6trhbZtkIjBVF8aDl77WCICQBi4EHAPsuXW/5GQPo3JR3QhTF5cBygPz8/KCbjxwuaeBsdYvD9Iju2FNQy7eHy3lg7tCgTtXOVDWjUSlC7p4f6qQURbVa5rjxNp41LDQd0aYTVTiMBREihmYmkBCtYsfZGq6ZEHGTYqcU1WnJTo4OS5jjp64aRbQ6si3QT1U2IYqEJVuWLS6/JVEU62wUvpWtQF+kpZsMQRBsNaZ1V7MFqHBTHnb2WMLYjnMRWdMZB0saeGX9KUqCnEj82dVHueG1H4N6DXue/uoIF/1tfciu16I3Ut2sd2v1oTOaOFbWSH1LcFM6VjYFPyG6K5QKgfx+KWFJEekPxbXakFvuWEmKjQqbh6unNOmMDEyPY2hmBC/vCIIwUBAEezdIq4b8EVACk2zKxlt+lgAH3ZSHnb2FdWQlRZOZ5L2VyvC2zdzgem0eLAndJq6VWLWSwtqWkHkBx0Qp2fv4Jdw4pZ/Leicrmpj7wgY2nawKansqGlrDYrljy+T+aZyubKaqSRfWdnhDKDNm2VNY08KfvzgU0bb6Uweksfb+GQwK8/KOu/nQr4Fn7Y79DCgQRbEC2Azcb1N2J3BYFMUyURSrXJX71+zAsKew1iv7fFuGWJT+kSDG1q9t1lNUq2V0iJV+39RYRJGQBbESBIHkWDVJsa43j62jyGBu2DXrjDTrTWEd6QNMHZBKfr8Uapr17itHAHqjmdpmfcgds6y06E28tfksR0q7XpL0UONuTf8tYJsgCEuBncCVwBXAbZbyh5Bs8r9DWgWdBdxgc7678rBhNJmZOTSDCV544tqSGB1FbkpMUD1zrbFfwqH0QRo9hSLRw8YT0qblvbMGO7TRt5IUE0WCRhVUs80mnZHBGfFOYwCFivF9U/jozvPD2gZvUKsUHH3qUgym8ISPyE6WZuuRHHjtule3MGNoOvfMGhzWdrhb09+DpKQvRdpgHQhcLori65byzcCFgBFpc/Z6URTftznfZXk4USkVPHnVKBb4YQo5LDOB0iA+ZG2euEEMp+wIW6UfCn44VslrG08TpXS/ARjsuPq9E6P59ncXdUrZGC70RrP7ShGCSqkIeThqKwnRUSREqyI2xLJWb2JXQS2mCPh3unWSEkXxY+BjF+XbkDoFn8rDRXWTjuRYNUo/LA1evGEC0VHBsxi4alw2/XvFkRwbWmeX9AQN10/qQ/9eodlwKqxtITcl1iNrIclWPzJf7EDz3rYCnvzyEDsfuyQk3tj+sOF4Jd8cLuOhy4aHra05ERxi+WSF1XIn/GEiItvGKYjc9/4efrpsq18yYtTKoJo1ZifHcOmozKDJd4YgCPzftWOYNrhXSK5XVKv1OJjcbRcO4JF5w4PWlpV7ivjJv7aELauRLX1SY2g1mNkd4jDevrDjbA0rthcS7WJ5LtjkpsTQrIvMEOTHLSlW7fM/h4PIHj4ECZNZZF9hvWsv178NhuaKzsfjMuCBE4AU6vfxzw5y5ZhsLh7hwobcA1n2NLYa+GxvCRcP7+2TdZG/iKJIbYshJC71hTUtHu+tTBng0K8vYJwob2JPQR2xYQyIZWVC3xSUCoEdZ2uYPiQ93M1xSXGtlsxEx1nPQsWym/L9mrkHk+MVjaiVCvIiIAR0j1T6pyqbaNIZXVvuOFLSdsfj1ErWHCwjPV7jWul7IMueA0X1PPbpQfqmxnZU+j50IE5xIevJYZ/y0c4i9j8xJ6izmWadEYNJbEvK7o6GVgO7ztYyJjeJtPjAW9hUNOroFa8Ji4ORPXEaFSOzE7uEZ25RGM01rUSqwgfIS4vjmgk5Ye0UrfRIpb+nQJouexNZswP/vQZUGlRKNf+KrUc4ogEhF5QaUKntfrpRTCV7QBUNSrVU13LukcIKBMydbfR96ECc4kJWTnIMjZZQuR7tKfjYGcVpVBx+ci5GD+P3n61q5uZ/72D5TROZMzLwS18VjToyEsNrrmnLpLxU3vnxHDqjKaKdj4prtUzunxrWNhwsruel70/y8Lxh9EsLT/5ZZ9wwuS83TO4b7mYAPVTp7y2sIykmyvfExK11YNSDSccYsRFTsw4ObgeTHow6MHvhMbp8hsPDvwJ+FQ38PcrSGajddyAf3Nje0bR1Ih07kw4/XTBMVUquUEFJ4RmS+6S3y1GowNHI34/OSBCEjpY7LjqQnDslB3GPrTS87IwqGlqd5ukNB/NGZ5GRoMFoEonUvVyzWUQQcBs7KdjojGbWHCrjp5P6RJTSN5lFTGbRpTlyKInQxyi4XDshl/x+qb4vW9z2fduvn246w1NfHmbH/Re3O/SYzVIHYNJJncNzLrLkXP+e1FFYOwzLz9fWHyUzVuDKUWltHQxGHez5r3NZVSfbr9nhpw7wzn562tfz2KQBVtiXCI47EVd8sthp53OkQseRSh1XTswjSh0j1XPRgaQay8mJaqCqsgIMmZIcVwG+vOyMBvdOYLCzRBwhWlqzlTWxX0pIcjb4g0IhsOkPs+gYWzH0WJeXIs1s82BxPde+uoU3F02KiL2ZHqn08/NSyc8LzFR0ZHYiwzITqG3Rtyt9hQIU0RDlwQbssMs7HWrWGXnmi6/57ZQhMNvOkcOV0r/bSYweUQSz0a5z0cE/xzoV1Tp/OY99vIsrR6Zx0YCkzp2IXSdFzWnn7Sr4sXN9y2xouOXDKuen2yK8MJrNSmCv5QPS7KPTTMb9bIbP7+s0G3qxj+W8HT90njW56kCaKtuvr1S77ois53h4vKZZz5mqZufKPwydkSPCGSMeJFNjlUKIOKV/vLwRo1mMmBlkj1P6hTUtnKtuYVL/FOdrpPpmJAdiByOXuI4hjqcOSGPNb6a7vmhchvMXyVF1jYo9j88JXJ5aQQBllPTxkOgJP6V//UTiB6RCPw86yIMfOS/7zf7OxyyzoXvf2UpFbSMf/GpCe6fwylTnsua/yL83nsCg03Lb+TmuOyKTHioOOZd1fE3H+qIf5n72szlFVOdlNtvfXfHVgx32hLYdqWFfmZbx88agUNkv27meGdFa3955eaKUfVim+/5oOe/+WMBfrxvTcXM9xJ2RUiGQlRwdcbb6JyqaUKsUEbPk1OOU/hf7S/jrmmPsffwS50p/138AEW75BvpO8f+i3j7g4DyJiZcdiEvcyLo7mMmbLbOh43VK+qRlQ5KHIYQn/ILJmQ2olAJ4YvP8hIsQFr8/3uHPQ0U13PfONp69aij5uXGdZzevzXQua95z7TMoR8tr9h2SK/a/336uaOYy4DKAr9zfbif+z2bz0FFHZN8hueKbx9pnUDazKvPxOuKPV5NwuhA0MZ7NjAytnndE1nM8OD4iK5FodxveIe6Mjpc3MjA9PmKsi3qc0t9TUMcAV16uRj1sfQn6XeCxwn/mqyMcKmngnVsD0EEAr288jVIhcPMF/TsX+tCBOMWNrHqtgeJarcuE5W340BmJokhRbQvnD/LO9t6j9vhAeZOBU3VmlHEpkODlOvrk29zXscVVZ/RQQfvvJiPF1XXMe34tj84dwMJxGR07EJMe3pzrXNacJZa6BvcdkclNcLftr0t17GZEFwMXq4FP3N51O3+xmDi7mxF5MjNa+2RbB7RskKVD2r3XuTWdqw7EbHa/NGd/jpvjJ8qbyM/z4HkKZGfkgh6l9EVRZE9BHdOHuPA0PfA/aCiGK5d6LhfYfrYGo8kcEDvc93cUkpcW51jph5BlP5xi+YbTHH3qUvf35cND2aw30TspmgH2VlRuOpCKhla+O1LBxSMyyEhws2/iRWdU0SCNwMMdYbMDShU5Gb2IS+rFD6VKFs50HX66E+ff6119V53RY5bguGZTh1nNr9/dhl6n5dXrR1k6EEvH8p8rncua/bilnr7jPpOjWZLRTc6KTS/4tzRny5Mplj0i287HptNQemFNt+4ZRKWa53IaSYmJh137OndAbddxs0wXQHqU0i+u01LVpGO8M6css0l6gDJHw6DZHssdlpmA3mjmbHWz37Gym3VGTlU2ccWYLL/kBIK+qbEYzSKl9a1BMceL16j4/v4ZnQvcdCCFtS08svIAWcmTyBjqRul70RlVNLpR+iFcWrNnUv9Utp6qDnlGM4colKCOBaRnYn9TIiOyc6H3SM9lXHi/+zq2uOqM/lTT1hFtPVHGq2sP8bcFw+gdp2jvTEyG9t/fuda5rBmPOJgJGRzPjlqaXbf5h/9DAM6z/r3Hu1sOFj1K6VszZY135vJ/9EuoPgHXven5WiMwLFNabjhS2ui30j9c2oAohj6csiOsir6w1n3+2lBijasf6MBrFY2tJLvKwBTCpTV77p01iHudheQNY2cEkJkYzYis0EaC7YSlIzJEJbChRME5Uxq903yw0JvxB+/qu+qMHq+lvK4B0aijd6wCwWQ/i7FdajPAip96314f6FFK/9JRmXxxzzSGZjpQzKIIG5+H1AEwYoFXcgdmxKFUCBwta/A7JO9BSzjlUGfLckSHEMsDAy///e0FfL6vhH/fPNkrx5X0eA1qpSLgcfWH9E5g3uiAigwYLgcTYeyMAFYsdmJtFYbOKLvNVr8FCK+HMAoFr20p4b8/nuPwk5fKG7nhIEqpYHSuE2V6ej2U7oUr/ymNGrxAo1KyMD8woYibWo30S4uld2Log6zZk5UUjVIhUBCkuPoHius5XNrgtaeiwmKaV1Qb2Hb94ry8gMoLNF/uL6FFZ2LhpD7hbopnhKEzsjpoldS52AcIYWd0vKKJQRmRY7kDPUjp641m/m/1Ua6ZkON4FL3peYjPhLG+JfZ65prADBHvnT2Ye2YF0VTSC1RKBc8vHNu2fBVoimq1Pjus5AYhmUpErJe74NM9JZysaIwopf/NoTL+8d0JXv9lftgDroEU7jw1Tu362QhhZ3SivJGpnkaGDWRn5IIeo/SPlDbw5uYz5OeldFb6RbvgzAa45Cn3O/Iu0OpNqJQCUX5a8ESS4rlqnO+ZxdxRWNvCUB/ji//1urHEBTBLkyiKjHniGxZPH8C99l7QEcLk/il8d6RcStweATNBgFOVzRwpbSAxOnJUyfTBvciMgO+nodVAaX0rgz1NnBLIzsgFkREBKARYI2s6DKe86XmITob8m32Wv+VkFSP+tIa9hXU+y9hbWMf8lzZxuCRykjsX1bbw7eHygMsVRZHiWq3PG8Q5yTEBzSjW0GqkUWcMW7o/T5hkCR2y/WzkhFourmuRchdHe+7tHWxeuH4890VAx32ivAmAIX4adwSanqP0C+vonaghyz4hSeUxyWpn8mLQ+P7P6Z8ehyjC0VLfFfa+wjr2F9WTHBs5L9Bne0u47e2dNAc4k1Sz3sSkvFSf8/+eqWrm798co6LBjQ23h1Q2SnIiykbfjlE5ScREKdkRQfH1i2vDH0c/UunfK45/Xj/O9xDuQaLHKP29hXWM65Pceelk0wugioEpt/slPzMxmqSYKI6UNfos42BxPWlx6s4dUxjpa2O2GUjiNSreuXWKz8tH5Q2tvPj9SY5bRlP+4tZGPwKIUiqY2C+FkvrAdHSBoLhOS06EBBKz8tneYib95TvqWtx4GAeZ1Dg1V43LCUqyH3+InIW4INKsM9KiN3W2z68rlDxwJ90Kcf7lgxUEgaGZCX6N9A8U1zMyJymi1vTbzTa1QdvQ9YWcDqZ5/lNpUfpuPXzDzBuL8iMqmcq4PskMjaDnAkCtVFDZqKO4ThvQJUBvWX+sgsyk6Ih6b6CHjPTjNCq2PzKbX02zC2uw5UXp53n3BOQ6wzMTOFbW6FN0zFaDiRMVTYzOiawHxKr0A222uXzDKWY9tx6DyezT+ZlJ0SiEwDlo5abEcMPkvmHJR+wNkaTwQdpQ7/RehRnrzCPQfhze8oeP97N8g4uQ42GiR4z0wUF2puYq2P02jPkpJAfGBG7e6CzyesVhMJvReGnr39BqYM6I3kzpH9zE396SHBtFgkYlOWgFkNOVzTS0Gn22dIpSKshKignYiz2xXyoTPQkhHWZEUeSeFXsYlZ3EnTOC4DHnZVsiaVZqJbvNVj98Sr9ea6C8QccQH63TgkmPGOk/8OE+/v7NsY4Ht/1LCuR0wa8Ddp0pA9K4+YL+Po3GMhKiefXGiRGRWccWQRD49y2TuP2iAQGVW1jb4ndSiZzkGMobA7O+3aQzBi5/QRARBIGSOi1rjwTeospbvj5Uzrgnv+Fkhe/7WMEgLU6NWqUIazKVE+XSdzLEU3PNENLtlb7ZLLLmYBk1zTabOq0NsH25lLUqfWhAr1dY08LpSu83F1sNAYoSGAQm9kslKymwm3VFfphrWnnr5km886vAhLNe/PZOFi7bGhBZwWZyXir7i+rD/swU12mpazGQGhdZG5WCIPCTibkMDuMo22pgMDjCzDWhByj901VNNOqMHe3zd70lZRS68HcBv95Nb2zjOftZhQdc++oW7lsRIWH47Dhe3shbm88ELAeqySxSUue7N66VOI0qYMsLlY06ekWYlYUzJuWlojeZ2eeHT0ggKK7VEhOlJCWCTIyt/OXq0SzMD5/n8vHyRmLVyog0Z+32Sn+3fWRNQytsfRn6XwQ5EwN+vWGZiRwt9W66qzOaOF7e2LYWGWn8eLqaP39xmMomNxmfPKTVYOKqcTnk+5nw+2BxPff/bx/lAbDVr2jUkZHYNZS+NSHHjjA7aRXVtpCTEhOR6/qAz0YCgeA3Fw/mg8XnoYigmDtWur3S31tYR0K0qj1Rx74V0FQO034blOsNy0rgTHUzWr3nU+/jZU0YTGJEhFN2RB/baJsBIE6j4rmfjGX28N5+yalrMfDx7iLOVLmJa+4GndFEvdZAehcZ6SfHqrlqXHbYZybFAZitBYvXNpxm6GOr0RnDswSWHKt2HtwxzHR7pZ+dFM1V47KlHtdkhM3/hOzxMGBGUK43LDMRUZSmd55ysEQKpxyxSj8lsGabeqM5IEtFgTLNa7PR7yIjfYB/Xj+e6yf3dV8xiFwyojdzR2aGtQ3OSI6NwixCWRgc2RpaDSxde8Knvb1Q0O2V/j2zBrNkgSUC5uFPofYMTPudV0lSvGGYJVb/0TLPnbQOFNeTGK2iT2pkjpqso7mC6sBYQ7y87iQj//S139Nvq+eyv7b6MVFKfnfJEMb18W+5KdTojWZa9IENj+ENv7l4CDeEueNxRjht9Y+VNfL8t8c5F6SQ5P7SrZV+q8HUboYnilLIhbTBMOyKoF2zb2osr/x8AjOHeh4OdfawDH5z8ZCIXRuNjlKSmRgdsFAMRbVaEqOj/I5GGh2lJD1B47dXblq8hvtmD3acXCdCqW3WM/qJr1mxvTAs19cZTV4tYYaado/t8Ch9ICJt9KGbK/23Np9l3JPf0KQzwsnvoPwATPuNd9nuvUShEJg3Osur0Lezh/fmlgjzarTnwzvO46mrRgVElpR+MTCzmoHpcfi7X1fTrKe8oTVg1kmhICVOTUaiJmzB17acrGb44/5FlQ0mVs/qcCj9E+WNxGtUZEeod3e3Vvp7CmpJjVMTr1FJqRATc2D0wqBf93RlE+9tK/BIidQ26zlW1ogxjJYGntAnNTZgYYeLa7VteW795f3F5/H3hWP9kvHvzWc475m1dAHfrA5Myktlx9masHRW1qxlkRQc0BaNSsmdMwY6DqUeZI6XS9myInXm3m2VviiK7Cmsk0w1C36Egi1w/r2gCn4Apg3HK3lk5YG2yI2u+PZwOXNf2EBhmOOEuGNvYR1LvjzstzWEwWSmtF5Lnwiy+qhs0pEap4molHaeMDkvlepmPaf9tF7yhaI6LWqlIqItnv5w6TBmeLHMGigKaloi0hPXSrdV+iX1rVQ26qSeftM/ICYVJvwiJNceliUFTTviQcTNA8X1xGtU9PPTOzXYnKpo4vVNZ/zeGDOYzNwzazAXDPIvqqmVLaeq+NlrP1LhRziGigYdGREcUtkZk/pbkqqEYYmnuFZLVnJ0RNqhW9EbzWGJv7PhwZn88YoRIb+up3Rbpb/X4pR1XnwZHF8DU+4AdVxIrm214DnmQWz9A8X1jMxOjOiXB6BvWmDMNmPVKn53yRCmeJo31A06g5ktp6r98iGoaNRFdBx9ZwzoFcfDlw1jUl7orY4i2UbfynPfHGPGc+tDvvylVAgRlUnMnm6r9AdlxHP3zIEMPPYaqONh8m0hu3ZyrJQI5agbpW80mTlS2uA4UXuEYbXV93cZqrZZT02zPmAvotU0zx+zzcrGrjnSFwSB2y8ayKAwxHe5YXJffjopMs01reQkx6A3mqluDl0ylW8OlfHwJwfCHhfJFd1W6Q/NTOCBSRqUhz+BiYsgNrRhc4dmJrhd3jlZ2YTOaI5YpyxbMhI0qFUKv71y39h0hkl/+Q5TgHZNraZ5/ij9B+YO5ZoJuQFpT6hp1hn5/mh5x4CCIWBhfh/mj80O6TW9xRrWJJS2+ptOVvHlvhI0qshVrZHbMj8wmMxsO12NcdNSEJRw3t0hb8Nfrh7Nh3ec57JOn5RY3lo0KWDr28FEoRDokxJDlQeb064orG0hKykalZ82+lbiNCpSYqP8Ms27dmIu5w2MrDwGnnK6splb/r2TjScqQ3bNFr2R05VN6I2RbXGWnSxZFoVyXf94eSODe0eu5Q50U6V/tLSRe5avQdj7Doy7ARJDPyLJSY5xu64Xp1Exc1hGl1lPXnXfhTz/03F+ySiq1bYtFQWKyf1TSYrxbQ21odXA3sK6sHq2+sPwrATi1MqQBl/bfa6OWX//gV3nakN2TV/ITZaes1Da6h8vb4pYpywr3VLp7yms5RbVahSiES74TVja0KQz8tzXx9h6qtppnY93FYU9PK43REf5b6dfWON/8hR7lt2Uzx8uHebTuXsK6ljw8mYOl/ie2zicqJQKJvRLCakFj9UDOtI3chNjVDwybxhTA2Q04I6qJh01zXoGZUSuuSZ0N6X/t8HwRBK/+Hosd6q+QBDN8OIE6XiI0agULNtwivXHKxyWm8wij316kJV7ikPcMt/ZfqaGu9/dTX2LwafzWw0mKhp1fidPCSTWYGtdZbbliCn9Uzle3kRtiNb1i2u1KAQiPp+wIAgsnj4wZIYSVU06+qbGRlwidHu6l9JvdqxgnR4PIlFKBYMyEpzG1j9d2YTWYOoSm7hWapr1rDpQ6pfZ5jPXjGbWsMA6zKw5WMbM59ZT7UO8f6t9f0ZCZCswV0zKk4wUQrXcUlSnJTMx2u/YSaGgorGVg8X1IbnWsMxENjw4k2mDI3uPLvL/a12Y4ZkJTm31D1gexK5grmmlb6p/tvrRUUpumNw34PesEOBMVbNPa7cVDToSNKqAhZgIB+P6JrP61xcGvDN1RlGtts1UNtL553cn+MWb28PdjIhCVvpBZFhWAmUNrQ6n3QeLG4iOUjAwPTQOY4HAGiTN12ibhTUtHCyuD3gCcn9s9SubuqZjli0alZLhWaFz8Lt75iDumjkoJNfyl+zkGGqa9SHZqP/1+3t4ZvWRoF/HX2SlH0SGZiYSr1E5HIEeKW1gRFZiwEwXQ0FCdBQpsVE+j/Tf2XaOa17ZEuBW0Ra8zRd77MUXDuDxKyPXZd5TDhbX87sP9nLdq1v8CknhCRcNSfcqdHg4sfpxlNQF9zsRRZENxytp0Pq23xVKvNI4giCMEwTBIAhCnuVvjSAIrwiCUCMIwnFBEC6zq++yvLszbVAvDjwxx+Fyxn9umcwrPw98jt5gMyI7EV/Hk9ZlgUCPSJNiokhw0rm6Y2yf5LAE5Qo0lY06PtlTzK5ztSz97kTQrtPYamDTiSrqu4Byg3YHrWDb6lc16altMTA4DN7R3qLytKIgCApgud05S4GFwD1AEvCRIAgTRFE85mF5YInLcLxpGxeel9pV1Ea1ShHx1g+OePfWqT6fWxQEc00r80Zn0S/Ne6ugrw+VMSIrMaIsirxFygUrOUqJwDvbCnhnWwEalYJjSwI7zjpS2siNb2zj7VsmM31IekBlBwPr0l+wlf6J8shOnGKLx0ofuBtoM4YWBCEXuBX4hSiK71qOTQDuBxa7Kw9M8+14IHgjHF9ZvuEUx8qaOsR833yyirVHKvjNJYNJjODATIGmqFbLnOzgbFw/e90Yr89p1hm5/b+7eOiyYdxx0cAgtCo0bHxwJku+OsKX+0owixAdpWDuyEwevXx4wK9ltdHvKhu5vRM0/PP6cUzoG9ygdDvPSX4SqXGR/z57tLwjCEIO8BfgYZvDFwEm4BObY18Csz0s7xGU1ev46kBph1gz649V8M62c8QEwNkp1Gw5VcU1r2ymtN67kVOzzkh1sz6oDj3eBnGz5juI5JjwnpCRGE2CRoX19nUGMwkaVVDMUItqpP+7da080lEpFVw1LifoM7ltFue497YVBPU6gcDTNf0Xga+AVTbHsoHToijavv2FQD9BEJQelPcIhmUloDWYOmx+HiiuZ3hWYpewc7ZHFGF3QR1nq7zbzFUpBf598yTmjc4KSrs+2FHA8MfXeLXWbHXMykjs2kofJMegn0zMpXeihhlD06n0wWfBE4rrtPSKVwfEOztUHCltCFpsoqGPrSbvoVVsPil53r+zrYC8h1Yx9LHVQbleIHC7vCMIwnykUftwwLa7jAbq7KprASWQ7EF5p/gEgiAsxrL007dvZIdt9ZThFu+8o6UN9O8Vh9kscqi4gfnjIjtCoTOstvqFNS1eBSnTqJRB3TCN10TRajBTXKv1OA5Pd3DMsrLspnxAmu0EM9hXcZ2WnADHTgo2y344xa6CWjY+OCvgsjc+OJM/f3GItUcqaDWag7q0FihcDjUFQYgHXgJ+L4qi/Q6pDmn5xharQXqMB+WdEEVxuSiK+aIo5qenR/4mkScM7h2PQoAjFietgpoWGnXGLuWJa0tWUjRKheC12ebhkga+PVwesJDK9uS22ep73q7uEILBHkEQMJlFGluDY13z6OXD+VMXM3HNTo6htK41KM9eRmI0ZQ06Wo1m1CoFOmPwltYChbv1hSXAcVEU33JQVoG0hGOLdbekxYPyHkF0lJKLhqQTZ/H4tGZp6kqeuLaolAqyk6O9Vvof7y7i3hW7CZb/kHVj0RuzzcvHZPHurVNI9jFCZyRiNJmZ8dw6nvs6OAZywzITg74pGmiyk2MwmsW2Tj6QiKLIsbJG0uLUfHrXBfx8Sr+gLa0FCnfLOwuQ1uDtu8gzwBdAX0EQskVRLLEcHw+0ArXAfjflPYa3bp7c9vvk/qnsePTikKdwCyTTBvXy2uqoqLaF3JTYoC09pMWpiY5SeOWglZEQHdEjMl9QKRWMzkli1YFS/njFiIA6/zW2Glh9sIxpg3q12b93BWwHBIE2k95XVE+Tzsgj80YzIjuRJQtGBVR+MHD3RMxDUtTWz+WW45cj2d4XAb8FEARBBdwGfC9KGm2vm/Ieh+1tR3KSBXc8c80YHp7n3ZplYY2WPkG03BEEgV+en8fYPsken/P1oTI2n6wKWpvCxZVjsqlq0vPj6cCGWz5V2cyDH+3vcmGorZZGwYir/8GOAmKilFw5NjgGCsHA5UhfFMXDtn8LglBn+fWwKIoFgiA8AKwQBGEIkAVMAC60nGt2Vd6T2HWuljve2cW/bpzIk18c4obJfbl+cvfYqPaUotoWJvYL7rLAw5d51xH949vj5KbEdonMZd4wc1gG8RoVn+8rDmjER+t+SVex0bfSLy2WD+84jyEB9pZtNZj4fG8JV4zJiuhE6Pb4NfcTRfF/wJVI1jgtwBxRFLd4Wt5TyEjQUNmo47sj5ewrqsfUxSc6m09Wcd4za51GELWnXmugodXYFrAtWIii6FVM+cpGXbcw17QnOkrJnBG9WXOwDJ0xcAm6rUtnXU3pa1RKJuWlkhQbWMUcHaXkk7su6DLB56x445GLKIpnoWPoFVEUv0Ky4Xd2jsvynkBuSgzxGhUf7iwC6LKWO1biNCpK61spqGlhaKb70VOcWsk3v50e9A3TZRtO83+rj3L4ybnEql0/2gaTmZoWfZd3zHLGrRcO4OoJOagUgVvTL67Tkhit6pJe5OuPVdCiNwXcT8ST5z/S6HreQV0QQRAYlplAVZMOlULoEvE5XOFtXH2VUsGQ3glkJAZ30zTLsknnyWZudZMeUewejlmOGJGdyIWD013Gf/KW4tquZ6Nv5Z0fz7F0beDCtBwvb+S+FXso9COhULiQlX6IsI4IopQKGoJkQx0qUmKjiFMrPX7gt5+p4b9bzwbNRt9Km62+Bxt2Vses7jrSByiobuGZ1UcCFkv+2evG8PLPxgdEVqjJSY4J6Ebu+9sLWX2wlNgumHxHVvohwhp/XGswBTX0bSgQBIE+qbEeK/2vDpTy1zXHgmajbyUn2fO4+sMyE/n+/ou88iruahTXaVn2w2nWHglMutBe8RoGpEd20m9nZCfH0NhqDMiAS2c08cmeIuaMyCStCw4avFrTl/EN29C3ENzQt6Fi3ugsj5cOimpbyEmJCbqZakaChiil4NGITq1SdFkF5imT+6eSkaDh830lXDnWv7Afja0G3tx0lstGZ3bJ5UmrX0FpXSuJmf7tSXxzqJy6FgM/ndQnEE0LOV1a6be2tlJZWUlraytGY/DTofnKxzf0o15roNVgwixKOV2jo5QkxURx5Ejkp1dzxCUWHeKs/VFRUWRkZJCYmEhhjZa+PsS69xaFQuCBuUMZ5UH45k0nqjhW3sgtF+R1aZ8JVygVAleMyeadH89RrzV4HJPIEeeqW/jHd8cZ3Du+Syr9dgctz4wPXPHBjkJykmOY1kVNfbus0q+vr6e8vJz09HQyMzNRqVQR/fIW17ZQ3axHIQiYRZG0OHWX3RSzYhZFBDo7momiiFarpbi4GFEUKaxt4fxBoVlGWTzds7j4aw6Vsmp/Kb+a1j/ILQov88dl8+bmM3x9qIyF+b6PTK2zp2CGxg4mI7IS2fSHmWT6aUwgiiKjc5O4eHhGyHISB5ouq/SrqqrIzc0lNrZrKE6jWSQtTkNqnJqaZj1Gs9n9SRFMU6uBM1UtDEiPI07T8TESBIHY2FhycnIoKCqmRW+iT4g6uMZWA0W1WoZnJbqsV9Gg63YhGBwxNjeJYZkJ1Hjhv+CINhv9LhR+wZboKGVbLmV/EASBP1w6zH3FCKbLKn29Xk9MTNd5APulxbX9nqPuOu12hkqpQEREbzIT56ROTEwMosnIgSfmhKxdb246yz++O86xJZeiUTm3rLAGvuvuCILAV/dd6PeotKhWS3SUgtQ4dYBaFnre316AQiH4POMxmUU2n6zigkG9AmoKG2q6tPVOJC/ndHfUlkBeeqPzGYv1/5MQHRUyN3Xr2m1pXavLepWNOjJ6gNIH2hR+k873fa+yBi05ycHfjA8mn+4t5oMdhT6fv+F4Jb94cztrj5QHsFWhp0srfZnwoVAIRCkVLpU+SPFJnl1zNOg2+lba4+o7t+ARRZHKJh3p3dQxyxF3vbuLW97a4fP5L90wgY/vPD+ALQo9OcmxfiVIf39HAb3i1UFNBhQKZKUv4zNqpQK9yb3SX7G9IGTT4faIis59CARBYN/jc7ini8VM8YfhmYlsP1vjs9JTKASSY7vu0g5ATnI05Q2tGNw8s46oaGxl7ZEKrp2Qi1rVtdVm1259D6GpqSkiTVJT4tRu4+kYzWLINnEBMpOiUQiuR/oAMWpll4qM6C9WO/0v95e4qdmZFr2Rhz/Zz65zgQ3VHGpyUmIwi1BW73rpzxGf7C7GaBZZ2EVt822RlX6Ec/bsWXr16sUnn3wS0uuWlZV1OlZdXY1e324FkhqnduuRaDKLITXzi1IqeO4nY7lslPPAWsfKGlny5WG/pvpdjbxecYzJTeKLfaVen1tcq2XF9kK3HWmkk50cgyC0h+Dwhg3HK5mcl8rAbuDQJyv9CCcvL48xY8awfPlyt3Vnz56NIAhtH6VSSXl5OTNmzOhw3Pazd+/eTnJ++OEHcnJyWLt2LUajsS35y+zZs7ntttva6plMJpq1OsxO1utFUZRG+qmhNau9ZkIuI7Kdm2weKW3g9U1n0BoCF3a4KzB/bDYHius5U9Xs1XlFXdxG38rUAWkce+oyJvZL9frc//5qCi910bhD9shKP4J47733HCrmHTt2sHbtWodl7733Xtv5arWaO+64g9raWu6//36mTp1K7969WbNmDSaTCVEUO3x0Oh1jxozp1I6XXnqJiy++mAMHDhAVFYVCoZDWwfft4+233267tkql4me/uJlmJwG9TGbJeSvUyqKotoUNxyudlltHej3FesfK/LHZPL9wrNf3XdRmo981fGKcEaVU+LQeL4oiSoUQ9CixoUJW+nZUNLSycNlWn6aA/qLRSC+jvXIWRZGmpiZ0Ol3b37W1Uprh6Oj2B1GtVqPRaEhOTubbb79lwYIFbXUUDuKqq9XqTsc3btzIRx99xOOPP84999yDwWDAbDYjiiKpqam88cYb7W1qaeXxZ19wasEjJVGP4edT+gXi6/GY97cXcvO/d2B0smFX0aAjOkpBvKbLuqn4REZiNNdMyO3kTOeO4lotUUqhW3SSz319jP9uPetx/Xqtgel/W8fXhzovd3ZVZKVvx9K1J9hxtiYskTCVSiX9+kkK8pVXXuHpp59Gp9MBsHbtWpKSkjhxQmqXSqXioosuIjGxfRnDakN9+vRp9u/fz9VXX01FRQVVVVXU1dV1+lRXV1NcXIzJJC1zNDc3s3jxYgBycnJQqVS88sorTJ06lWnTptHS0sLzzz/PtGnTmDx5Ml99+TlqjcatBU+oHVlyU2IwmUXKGhx33JVNkjduV7Y595V6rYHXN57mZEWTx+do9Ub6psZ22bADtvxwvJJvvYg6+vneYgprtF3WE9kR3XKo89NlWzsdu2JMFjedl4dWb2LRW9s7le86V4vRZm3aGglTEGByXio3Tu3HlWOzKanT8tsP9nY6/7YLB3DxiN6cqmzyebNnwYIFbaPzwYMH8+CDD/L666+zdOlSFAoFOp2Ovn2l3Lrx8fGsX7/eoZyPPvqIUaNGMWjQIMaNG8e+fftcXvfMmTP07duXRYsW0dLS0dSxqqqKrKwsPv300w7HZ8yYQWNjo2S26WSkX9usp6ZZj9kshlRhtAXXqtU6dL1v1hm7xajVF4wmM8+sPkp1s97jcAJ/vmpU275OVycnOYaTlZ51eKIosmJ7ISOzExnVxbPd2SKP9C08fNkw5o/LJtqy5qcQIC1Ozbjc5JC2Y9myZdx8882MGDGCnTt3cscdd6DX62lpaSEzM7NtCcgVI0aMoLCwkNLSUrZu3dq2GWt9cdetW4coiphMJpqamujXrx9nz55lw4YNvP766x1kKRQKVq1aRXJycofPpk2bUCikNVJnSr9Zb0RnNId8hGgdlTmzNnn9l5P44PbzQtmkiCEtXsO0Qb34Yl+JV4q8u8yKspNjKK7VenTvB4sbOFzawPXdwEzTlm450nf1QseolU7LH115AJ3JjEYlOR1dNiqTJVeP7lAnOznGpXx/TbqmTZvGl19+yfDhw1m6dCkPPvggAC+88AIDB3oWQfKKK64gPz+fJUuW8PLLL3P69GnS0tJISmofrVhDUvfpIz3QAwYM4MiRIw7X/i+//HKHI32AXvFqzE5eIL3RHJYYJdltDlrOTQy7cuwUf7lybDa//3AfewrrmNA3xWVdndHE4rd38Yvz+jF7eO8QtTB45KTEoDWYqGsxkOImjtD7O6ScF/PH5YSodaFBHunbUNWk4+dT+rHyrgv4+ZR+VDbpQnr95uZmhg0bxhdffMFLL71Ebm5uW9nZs2cZMGBAh/p6vZ7mZsfmd1dffTWfffYZAIsXL2bhwoUdyq+44gp+/etfdziWmurYlG3NmjVkZmZ2+GzZsgWQ4uokxTh+efQmM6owKNfoKCVv3zKZ6ybmdiprNZi4+73dbDzh3LqnuzNnZG/UKgWf73XvqFVa18oPxyv9jtIZKeSmxJCRoKGmxf39zB+bzWNXjPArD0Ek0i1H+r6y7Kb8tt+XLBgV8uvHx7ufJbz99tsd/rYuzVgxGAwcOHCA77//nsTERI4dO8b333/P2rVrO5z3yCOPsGDBArZt28aUKVNcXvPSSy91OtI3m0W0BhMalQKVsn0MIYoiBpMYthH19CHpDo9XNupYtb+UiwY7Lu8JJEZHcfHwDI8UuXW2lNPFbfStzB2ZydyRmR7VnTIgjSkDul86TVnpRxBFRUVER0ejVHYMCfzWW2/x9NNP8+OPP5KYmEhUVBSiKKLX6zHbxeUXBIF58+YRFRXFihUrWLJkCTNmzGDmzJkd6l111VWcd955PPTQQ6xbt85lu6xr+rY0NTWxaNEidEYTpyqb6JcaS5JNbBaTWZSWyZThUfoHi+s5Xt7INRM6jvYrGqXZW08KtuaIF2+Y4FGHbI2jH8pQGpHAm5vOMH1IOoMyur4Hrj3y8k4EkZOTQ1paWttmaUJCAm+99RaPPvoob731Fi+++CLnnXcea9asITk5md69e5OV1R5uwGg0olKp+PHHHzlx4gRRUVG8++67PPbYYw6v98gjj7B+/Xq+++47l+26/PLLO5l7Tps2DaDN2UVnZ7apUioY0juBWHV4xhVf7C/hoY8PdPIWruyhjln2WBV+qxuv5KI6LQpBimnUXbhvxR5e33jaafmpyiae/PIw3x7u2iGUnSGP9COQs2fP8r///Y/XX3+dxsZGVq5cydy5c5k6dSrx8fHcdttt/PWvf+XZZ5/lkksuaTvPYDAAUucB0mh88eLFzJo1i1WrVrWZbqrV0oj88ssvZ8WKFVx44YWdZNTX11NWVkZraysGg4G6uroObTQajTQ3N1NeVobZoMJgjKwIjLkpsehNZiqbdPS28aSstI70e7jSB3h1/Sle33iarQ/PduqpqlEpGNsnmShl9xkfHiypd5m57n87ClEqBK6d2L02cK10n/9kN2DlypX079+f/v37s2zZMn71q19x9OhR5s6dC0CvXr1YsmQJx44do3///syZM4d//OMfbedbFbaV6dOn869//QuQHLZee+01brjhBsaNGwdIS0HXX399BzNQqzPYxo0bycvL45VXXmHz5s3k5eV1+Bw8eJCHH36YgQMHUnD2JDo7s82qJh2nKpsIl3l3rhOzTZNZpFe8hrQ4WekPzYynulnPppPON7XvnjmIlXddEMJWBZ+c5BiKnSTZ0RvNfLy7iNnDMrptOk15pB9BLFiwgKKiIi644AImTJjgtF52djYff/wx3377bdsyC8D333/v9Jx7772Xe++9120bcnNz22yY77nnHo/afa66udMygVZvQm80Ey7z7jYHrTotE/u1myUuuqA/iy7o3snQPWXaoHSSYqL4fG8Js4Z1fXNMT8lJjuGIE6/c74+WU9Wk5/rJ3cs23xZ5pB9BCILAvffe61Lh23LJJZdERJ7gjARNJ89XvcncllIxHLQ7aDlPptLTUasUXDYqk28Pl6PVd17bN5lFLvvnRj7aVRSG1gWP7OQYqpp0DvczSupaGdArjund2LpLVvoyfhOjVnUK4mUwmsOaYShOo+Lr30znl+fldTj+wIf7eHndyfA0KgKZPzabZr2J7492HvmWN7RypLTBbUrMrsbgjHjG9kmmodXQqeyWaf359ncXdTA/7m7IyzsyfmMym2lsNRKrVqJWKTGLIgaTOeybf0MzEzod++F4JYpuElIgEEwZkMYfrxjRYQnMSnez0bdy2egsLhvdOclOdZOOtHhNt/fW7r7dmUzIMJpFCmpaaNJJ02WzWSQ+OooYtdLNmcFlw/FKXtvQbppnMotUN+tlyx0blAqBX03r79Aks7gtjn73UvqOMJlF5r+0mT9+ejDcTQk6stKX8ZsopQIB2pYBVEoF/XvFhd19/YfjlTz/7fG2jemaZj0ms0hGD3fMssdkFvlsbzFbTlZ1OG7dD+luSt9sFpn/0qYOA4JNJ6sortMyZYD3WbW6GrLSl/EbhSAQpVS4jasfanKSpeBa1nAD1sQ46W7y+vY0FAL8dc0xlm3o6LCUnqBh1rCMsM/YAo1CIVDe0Mrx8sa2Y//bUUhKbBSXjOj+Vkyy0u8CtLS0dLLBB8lByhr//syZM20JVjzh888/59VXX3VY5klSdHtsQyyXN7RyrKwh7DHYbc02QRrRjshK7HZr1P4iCAJXjs1m08mqDvF4fjqpL28umhTGlgWP7OSYtueiuknHN4fLuGZCLhpV9+rgHCEr/QhiwoQJZGdnd3CCevPNNxk9ejRJSUkkJycjCAIJCQltoRqsOW6XL1/O+PHj+e9//8uRI0faErLYfqyOVwA7duxg9erVndrgaVJ0s9ncoSOyTaaiM5oRxfDHYLePqz8mN5mvfn0hY0KcI6ErMH9sNiazyFcHSsPdlJCQkxxDiUXpr9xTjMEk8tNuFjffGbLSjyBWrVrFnj172LNnDzfccAMxMTHMnTuXU6dO0dLSQnFxMQC7d++mrq6OpqYmTp6UzA+feeYZnnrqKV566SVqamrYsmULixYtYtGiRcybN4/PPvuMqKgoHnvsMQ4fPoxSqXQYO9/TpOhKpbJDJ5CRqGFQRhwgre2H01zTijVIWGl96PMddzWGZyUwMD2OL/ZJ4ZZFUWTSX75j2Q+nwtyy4JCTHENJfStms8jPpvRl+U0TGdK7s7VXdyT8b6ZMG1lZWTz99NO89dZbvPnmm6xevbotjg5ISzhxcXFOk6n89re/ZcuWLcTFxREbG8tFF12EQqFo89pVKBT85z//cbh8A94lRdfr9SxfvrztXLVKMtcEIsJcEyAxRsW+x+fwq2mSB+4L3x3nxte3hblVkYkgCMwfm0Ndi4FWg4nqZj2VjTo0EdB5B4OxfZKZPSyDVqOJWLWKOR6GW+4OdM//qC/8bTA8kdT587fBIW3GRRddxKOPPsovf/lL8vLyAPjyyy/blLdKpSI1NZXk5GQUCgVffvkloijy7LPPUltb2yEsc01NDVdddRVVVe1WGc5G+N4mRf/ss8/aAreBlHu1srGVFr0RgykyRvqCIJAU225BdLy8kdJ659m0ejp3zRzImt9cSHSUsm1JLKebhlSeNzqLV2+cyGsbzvDutnPhbk5ICf+bGSk0O47F4fR4gBFFEa1Wy9VXX83nn3/OyJEj28o0Gg19+/alrq6O3bt3s2vXLurq6ujbty/R0dE0NDTw2WefkZ+fz9GjR9vOs4ZosI+Fb4/ZbHaZFH3Tpk1otVoOHjzIpk2biI2NpaGhoWP7kZZRGluNpMapiYsQi4+Ve4p4do30nVQ26mQbfRdEKaWlPL3R3CNs9E9XNvHC2uNsP1MT7qaElO7nkbv6ISg7EFiZb13uXf3M0XDZ/3l1SkFBAYMHDyYmJobm5mbUajV33303BoOhLe0hwN/+9jdqamr44IMPAFCpVCQlJbF+/Xquu+46NmzYwOTJk4H2jVR38XmsSdHfeecd5syZ03bcNim6LdYEKraoFAIKQcBkFjvF4QknewrqWLmnmD9cOoyKRh1j5U1cl6w7WsF9K/Ywf1w20P28ca0064zM+vsPAOjc5BTobnQ/pd9F6devX5tJ5IwZM1i0aBEmk4m33nqL6Oh2b0mtVutwTV+tVvPpp5+iUCjYs2cPAA0NDajVaqKiXDtJ+ZoU3RZBEFCrFLQaTIiiGHbLHSu5KTE0thqp1xqoaJBH+u4YkplAo87IwZIGrpuYG3YHu2Aw9LHVHUKBrzlUTt5Dq9CoFBxbclkYWxYaup/S93KE3cYTSc7Lbl7lm0w/SU9P7xTeuL6+nvz8fIf1Kysrueuuu7j99tsBKC0tJTs7u63c1mTTntTU1E6JUqA9KbotNTU1nUb6IJltNrQaOFTSwMjsxIhQ/DnJ0qzjbFUz5w9MY1ROYphbFNnkJMcwKS+F6iYdBTUtVDS2dru48hsfnMmSr46w+kApBpNIdJSCuSMzefTy4eFuWkiQ1/QjmPnz53P99ddjMrVPP48ePdphpG9rK//aa6+xbds20tPTmTNnDtu2baNPnz7Ex8fzy1/+ElEUMZlMDh29nHHppZdSVlbW4XP++ec7rGvdvI1SChGh8EEa6YOUG/eNRZO4enyumzNkrhybzemqFnacqWHpd547/HUVMhKjSdCoMFryOOuMZhI0qm7XuTmj+430fSUuw/GmbVxG6NtiR2urZGd+6tQpTpw4wcSJEwEp4YnVWsdkMrFs2TJ+//vfM378eJYvX8706dOZOXMm27dv59///jcAer3e5YjfHk/X9AF6J2pobDViNJsjxmwzJyWGWLWSJp3nHV1PxnbpQwTe2VbAO9sKut3SR1WTjp9P6cfPJvflve0FbbmTewKy0rfyQOSOaC6//HLmzZvHtddey6WXXsq2bds4d+4cmzZtaqvzxhtvoNPp2swut27dysaNG1m8eDFXXXUVb775JjfddBPl5VKy5507d3p8bU/W9AGUCgUGkxmzKFLR0BoR5n5pcWoO/Xkuqw+WMeXp7/jw9vPpmxb+dkUq1qWPNQfK0JvM3XbpY9lN7UukSxaMCmNLQo+s9COI0tJSqqqqKCoqIja2XTE1NjZy9913s2PHDn788Ueqq6t5/PHH+frrr/nPf/6DwWDgscce46677iI2Nha9Xs/vfvc7rrvuOm688Uaqq6u56667mDZtGv37S45KjpZ4vE2KXlJSQlpaGhqNhoPF9ZhtYu1UN+upbtZTUaclnOrCusxUWt9KeYOOxBj5kXeFdenDYDb3yKWPnkD4598ybezcuZPZs2czePBgLrnkEgA++OADhgwZQnl5OTt37iQnJ4cxY8awbds2EhMTmTZtGmVlZaSnp3PHHXcAcPvtt3P06FGee+45AO677z5Gjx7NqlWrMBqNXHbZZSxdurTDJi/4lhTd6hcwNDOB5Bg1ApKSVQgCybFqMhPDryyWbzjFU18eRq1UdEtrlEBjXfpYedcF/HxKPyqbPF8OlIl8BE8jIQqCsAC4xvLnh6IofmE5rgH+AVwPVAG/FkVxtc15LsudkZ+fL7pagjhy5AjDh3evKacjtFotW7ZsYfbs2Z3KRFHk1KlTDBo0CJPJ1La+v3//fqqrq5k5c2Zb3cbGRhISpNgiy5YtQ6PRsHDhwg4zCn8prm2hulmPQhAwiyJpcWoays6F/f/08CcHWLG9gJzkGDY/NCusbZGRCTaCIOwSRdGxiR8eLu8IgnAX8BfgdSAZ+EwQhEWiKL4NLAUWAvcAScBHgiBMEEXxmOV0d+UyLoiJiXGo8EFauhg0aBBAh/AL1sibtlgVPtBm0hlojGaRtDgNqXFqapr1GM2REV/fasETr5GXdmRk3L4FgiAkAP8HXCuK4reWYyJwmyAI3wO3Ar8QRfFdS9kE4H5gsSAIua7Kg3A/MmGkX1pc2+85aknRHglNFAuXWJV+VbOuW9qdy8h4gydr+hrgHqvCt1BqOX4RYAI+sSn7ErAOTd2Vy8gEHWv8mJomfbe0O5eR8Qa3I31RFKuAt61/C4LQC7gReBPIBk6LomgburAQ6CcIgtJduSiKfgW9iCR3f5nOhDtzFvQcu3MZGU/xynpHEIRXgGPAfuCvQDRQZ1dNCyiR1v7dldvLXywIwk5BEHZWVla6bItarUarlcPkRjJardZt3J9gs/HBmcwfl010lPSoR0cpuGpcNhv/MNPNmTIy3RNvTTY3A7uAmcBUQIe0fGOLNclmjAflHRBFcbkoivmiKOanp6e7bEivXr0oKiqipqYGg8EQEaNKGQlRFNsyfWVkhNej2Wp3rjPKducyMuClc5ZlM/ZdQRDeAV4EXkBawrElxfKzBahwU+4zSUlJaDQaKisrqa6uxmg0+iNOJsBERUXRu3dvEhPDH+CsJ7vcy8jY44n1jhroLYpioc3hVcB1SMs8fQVByBZFscRSNh5oBWo9KPeL6Oho+vTpGcmMZXynJ7vcy8jY48nyzkzggCAIqTbHBgPngL1AEfBbAEEQVMBtwPeitN7irlxGRkZGJoR4srzzPZLiXiUIwp+BTOD3wMOiKJoFQXgAWCEIwhAgC5gAXAjgrlxGRkZGJrS4HemLomgALgcqgf8BfwQeEUXxZUv5/4ArkaxxWoA5oihusTnfZbmMjIyMTOjwaCNXFMVzwHwX5V8BX/laLiMjIyMTGuQomzIyMjI9CFnpy8jIyPQgPA6tHGoEQahEshDylV5IoZwDgSxLliXLkmV1FVn9RFF06t0asUrfXwRB2OkqprQsS5Yly5Jl9QRZ9sjLOzIyMjI9CFnpy8jIyPQgurPSXy7LkmXJsmRZsqyOdNs1fRkZGRmZznTnkb6MjIyMjB2y0peRkZHpSYiiGPEf4GZgvd2x/sDHQCNQBvwZUNiUjwe2AM1IKRqfQurkbgbW251fCRx2IcvhtWxkaYBXgBrgNLDJV1l29/yjm3t0167ewEqgCSmcdaE/7QIWIKXO/AQpoY5f92ipOw8wW9royz1ej5QJ0fZj9PO7jwbOAkd9vMejDtpk/dT7cI+2z3IJcMSLdj1n8wxogc+BT23K11t+FgA32/1vBIt8a/mv7WStAJ62Kf/MT1lxSPmzT9nV9VbWh8B/kbL26YGTSLk9fG3XXOB14D1LmV/3aKmbjfRO1ngoy758Kp2frZMe6dNwK3QPFP4UpEBtti9lHHAG+BqYAdwCNADPWcpjLeUvAdOBO5FemH9YZG2wOf9SJKVvBD5wIMvZtd6xtgtYhpQf4BYkhwojcJMvsuzuWeviHj1p10akTuh+S/vMSArBF1l3WWS8YCkTkV54n+7RUjceKbuaiKS8fWnX/yEpm2lAMbAVKYezP+1aYvkffufjd78JyLd8rO3aY/n+n/BSlvVZfQmYQ/vz9ZKH7TIiKb97kQY+IpLT4wzgI8vf3wE3IHUEM22+h4ct/5+7LeVGpE7nXuBxpKx4Rkv5CousF3yUZQD+A5QjPWe2db2VZbZ8dw8i5f4QkXJ7+CLLaGnbS5b/qwh86cc9LrXUPejl92Vffgewg/bnLB8Y1eWVPjALaWS0m44v5a2W44k2x36HpCTVSKOFekBpU77S8k/YjZTntx5ItJH1sM35trIcXesVyz9sD5KSMQE/t6n7H2C5g3a5k7Xe5p7PWdrr7B7dydqL9PJk2dR9DSjwQdZGJOVyiU3dt4CNvtyjzfH3aB+l5Pn4fa1BerndPRMetQtppKxHUhp+ybJ7Vl8DPvPx/1iPlFfaWvdfwDp3spDegxa7ch3SjE9j+Z9+ZVP+GPCN5Vxr+aM2srQ235PG8j3V2dT91uZ8j2VZjr2P9B7tQlLaj9qUedOuQRY5pTZ11yApXI2XshSW77Papu43QKGP9/gK0qj9atqf+xnuZNmXW/7+F/CKL3o10tf0LwR+gTQltWUCsFMUxQabYyeQpuW9LB8F0hTJSl+k3v5zpBGm9fwJwE6kntd6vq0sR9dKsvz83vK7CWnJwyprJdI/3b5d7mTZ3nMF0OjiHt3JagTOE0Wx1KZdZUCUD7IUwD2iKH5rI6sA6eH05R4RBCEHKfvacTrirazxSIrC3TPhUbuAfyK96FsDIAtL3f1Is48/+yArivZn2frdK5GUtztZO5E6BdvyEqTouvlAAtJSj7X8S2C6JdmRtfwDG1k3AFNtyqOQFL+17lKb872RhaW8CVhtuVdrXbyUVYLUwWJT9wfLd6bwUla05Vzbe1xL+3Pv7T1WI3Wur9IZV7Lsy6H9ufeaSFf6T4mi+JmD4yYgze7YcKTevBJpvVkJPC4IQoIgCDOBIUhTKZB6WOv5Vlm259v+7uha+y3lzUgPwGlRFLU2dQuBfoIgKL2UZXvPZtoVtKN7dCfLJIriUbt7vARpn8BbWQZRFN+2kZWBpMSs/xtv7xGkHMtnLfV9vUe1pS1/QRrFThcE4S82L4ZX7RIE4VKk3A9aYIIgCH8UBCHOj3vEUncwsE0Uxd0+yGrA8iwjva/9kJbC3nMnSxRF62zWtjwD6Rmwrikn25RbZwC5NuWnbGRttCsXkd41a90tNuXeyMLS7u+RcmiL1roWPJYlimILMM7uHscC+yzvqLey7L+vq2h/7r29x0uQBmP1lrq2OJVlX27RK6OBOwVBaBAEoUoQhFcEQYjHE3yZHoT6g7QOajtNWoj0YPzK8vdwpJHxFzZ1fkbHTY5XbWQdsp5vI6sB+MJelqtrWWSdAX60q/uY5edUL2XZ3uOHru7RS1nWuiJS+kt/ZH1tqbsbqVPyWhZSboZqJEXdtrzjg6wDSIrsGeBJS10D8KiP7dpD+5q3aJG9F0mJ+Pp9WTea/+Xr/5HOz/I6b2VZ/l5iKf8r0p5TpV15rKV8oqW8zO49tC1/lvbn6SakWaRtuTey5tjIWo40YHFW1xtZ1ns0YtkI9UPWx5bnYSuQ4IesVuAPlu/LdnnHG1mjLb8vR+pIbkNayn3NI30aboXuo9JXWf4JIu0blCJwuaW8F9Lm2YdIL8LjSGubv7fI+sHufOvL1OhAltNrWWSdAjY7qCv6IMv2Hp+0PLDO7tEbWfFIU+dAyLoJabPN1+9rIxZrBEvd1bRbtXgraxOWzSu7uiYfZO2i/UUK5P/xOto3q3357jfT/ixfT/uApcVLWTFIS0FWZWN9HuzPF5E2nxcCxXbvobV8FtKaeatNZ1Nsd743sk4AX1mOvwwYndT1VtZCy3e/D1D5KesZpBmgFljoo6w6pCU+6/dlq/Q9lTUNacN+vF35zZZnIq5bKn2b4+cDiyxf4C6b4w8gjdoEm2N3IimpJ2nfqLE9/6gjWa6uZWnXUeCMXd1bLP+gEi9lrbe/Z2f36KWsF5FMw67yV5ZN3TVIyttbWUXAdzb18izf1e8D1K7f+vjdH7acN8yubgmSgvT1u/8AadPb1//jKTo/y39DUj7eyLI+A7mW8v+z3K/te5NuOTYBybrHQEdjCGv5+0gzYwPS0pO1bm+b872RVYO0nIJNuxzV9VaWdUN8QgBkWev+BcmCSuWlrGak/Su1TV1bpe+prAn2/29L+TBL+Ri3+tRTxRvOD05ecEtZvuVmr7A5tgxYYVdvmqXec3R8KdvOdyTL1bUs7dqJNLLMtqm3GGka562sTkrf2T16Kgu4xlJ2rR+yfgD62NW7gfYRrDeyrN+Ls483sn4Ehjholy+yDliOaezqHvVBlvW7j0YaaMz347svwfmz7JEs+2fAUj7bcux+m2PWJYhMpNmyGZjqoFxE6lTMSEuY1rr32ZzvjSzbdt1jOeaorseykMye9S7qeiLrJ7RblVnr3opv93jQ5nf7z3ovZGUirfmPtft/n2cpn+xWn7qrEAkfXCv9/2FZXrE5tgTpZVXaHWu1/Fzv6HxHslxdi/bR+Dngb5ZjKiT72RJfZDn62492jUNSOi+6+r48kLUPaWqaanP8caTRnreytlvaZf3MszysG4AdXsoqAlba1fsCaVTlbbu2WF60yTbHM5DWg4/78n9EWlbRAzF+fPfncPwsm4At7mQ5egZsyuuBD22OrQQO2fy90a78e8t1X7Qvt/xebD3fW1k2xxchvaf+tKsP0vv3uYO6HssChlr+/4Nt6m5Heu7VXsoaRMfnfjfSc3+rpcydrLZyJIus3XSc/b2ItLwT6+yZaKvrrkIkfHA+lR9t+VIvtDt+PtILvAFpLe4jyz/vX3R8KdvOdybL1bVoV64LLWWfWR4Kk6+yHMj2tV0/IPkjFFi+j3zgp5a6M32QdRBpE+tSpBezyfId+3yPlmN5tK/DeytrJ5JSXYpk5mrd/PapXcC7lu/saqTOaJ9F3mxf7hH4O5aOzI//424cP8ue3OMMB8+A7XNwp+X7+w7JFFEErreRc4FN+feW8iobWb9CWoZYa/lfiEgzJl9k5dvILPWzXTuRns9Lbdq4zfI/9lbWViQrs+uARyzlZ/28R7Wlroi0l+ROVodyIAepw34XaY/tX5bypzzSp8FU1oH64FzpfwJ86eScKy0PoB5pU+ttJNvXNlm257uS5azcTtY8JOW4HqmX9lmWnSLytV3Wl9DRJ8/bdiGZCn6ONMo5hWTV4tc9Wo7lWdq01hdZwLVIa6WtSApys6/tQjKJ+yuSeVwz0kb6Oj+eiZ20jz59fr7o/CyfA1a7k4W0PuzsGVhrqTsFaX9mC/BTB7Ks5ftcyPrBcv6jAZD1eySlGoh22X8O+iGrHmlm+RqSBZs/7cqz1BWRBonuZHUqBy5CegdbkcJyLHL2PNh/5NDKMjIyMj2ISHfOkpGRkZEJILLSl5GRkelByEpfRkZGpgchK30ZGRmZHoSs9GVkZGR6ELLSl5GRkelByEpfRkZGpgfx/7cSKImcO92WAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "'\\nGM(2,1)\\nimport numpy as np\\nimport sympy as sp\\nimport matplotlib.pyplot as plt\\nx0 = np.array([412.0,559.2,380.8,542.4,553.0,310.0,561.0,300.0,\\n            632.0,540.0,406.2,313.8,576.0,587.6,318.5,400.3,590.6])\\nn=len(x0)\\nlamda=x0[:-1]/x0[1:]  #计算级比\\nrang=[lamda.min(), lamda.max()]  #计算级比的范围\\ntheta=[np.exp(-2/(n+1)),np.exp(2/(n+1))] #计算级比容许范围\\nprint(rang)\\nprint(theta)\\nx1=np.cumsum(x0)  #累加运算\\nz=0.5*(x1[1:]+x1[:-1])\\nB=np.vstack([-x0[1:],-z,np.ones(n-1)]).T\\nu=np.linalg.pinv(B)@np.diff(x0)  #最小二乘法拟合参数\\nprint(\"参数u：\",u)\\nsp.var(\\'t\\');\\nsp.var(\\'x\\',cls=sp.Function)  #定义符号变量和函数\\neq = x(t).diff(t,2)+u[0]*x(t).diff(t)+u[1]*x(t)-u[2]\\ns=sp.dsolve(eq,ics={x(0):x0[0],x(5):x1[-1]})  #求微分方程符号解\\nxt=s.args[1]  #提取解的符号表达式\\nprint(\\'xt=\\',xt)\\nfxt=sp.lambdify(t,xt,\\'numpy\\')  #转换为匿名函数\\nyuce1=fxt(np.arange(n))  #求预测值\\nprint(yuce1)\\nyuce=np.hstack([x0[0],np.diff(yuce1)])  #还原数据\\nplt.plot(x0)\\nplt.plot(yuce)\\nepsilon=x0-yuce[:n]  #计算已知数据预测的残差\\ndelta=abs(epsilon/x0)  #计算相对误差\\nprint(\\'相对误差：\\',np.round(delta*100,2))  #显示相对误差\\n'"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#15.2\n",
    "\n",
    "'''\n",
    "GM(1.1)的预测效果非常不好,误差较大\n",
    "GM(2,1)的预测效果极差，误差非常大\n",
    "知乎上有朋友做二次指数平滑，应该比较合适\n",
    "'''\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "from matplotlib.pyplot import plot,show,rc,legend,xticks\n",
    "rc('font',size=16); rc('font',family='SimHei')\n",
    "x0=np.array([412.0,559.2,380.8,542.4,553.0,310.0,561.0,300.0,\n",
    "            632.0,540.0,406.2,313.8,576.0,587.6,318.5,400.3,590.6])\n",
    "\n",
    "print('*'*10,'级比检验')\n",
    "n=len(x0); jibi=x0[:-1]/x0[1:]  #求级比\n",
    "print(jibi)\n",
    "bd1=[jibi.min(),jibi.max()]    #求级比范围\n",
    "bd2=[np.exp(-2/(n+1)),np.exp(2/(n+1))]   #q求级比的容许范围\n",
    "print(bd1,bd2)\n",
    "\n",
    "print('*'*10,'建模')\n",
    "x1=np.cumsum(x0)  #求累加序列\n",
    "z=(x1[:-1]+x1[1:])/2.0#累加序列的均值生成序列\n",
    "print(z)\n",
    "B=np.vstack([-z,np.ones(n-1)]).T\n",
    "print(B)#获得最终的B值\n",
    "u=np.linalg.pinv(B)@x0[1:] #最小二乘法拟合参数。求的最小值u的估计值（即使Y-Bu=0）\n",
    "print(u)\n",
    "\n",
    "print('*'*10,'解方程')\n",
    "sp.var('t'); sp.var('x',cls=sp.Function)  #定义符号变量和函数，x(t)\n",
    "eq=x(t).diff(t)+u[0]*x(t)-u[1]  #定义符号微分方程\n",
    "xt=sp.dsolve(eq,ics={x(0):x0[0]})  #求解符号微分方程\n",
    "print(xt)\n",
    "xt=xt.args[1]  #提取方程中的符号解\n",
    "xt=sp.lambdify(t,xt,'numpy')  #转换为匿名函数\n",
    "\n",
    "print('*'*10,'求预测值')\n",
    "t=np.arange(n+1)\n",
    "xt1=xt(t)  #求模型的预测值\n",
    "x0_pred=np.hstack([x0[0],np.diff(xt1)]) #还原数据（预测数据）\n",
    "x2007=x0_pred[-1]  #提取2007年预测数据\n",
    "cha = x0-x0_pred[:-1]; delta = np.abs(cha/x0)*100\n",
    "print('1989~2006的预测值：',x0_pred)\n",
    "print('\\n-------------------\\n','相对误差',delta)\n",
    "t0=np.arange(1989,2006); plot(t0,x0,'*--')\n",
    "plot(t0,x0_pred[:-1],'s-'); legend(('实际值','预测值'));\n",
    "xticks(np.arange(1989,2006))\n",
    "show()\n",
    "'''\n",
    "GM(2,1)\n",
    "import numpy as np\n",
    "import sympy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "x0 = np.array([412.0,559.2,380.8,542.4,553.0,310.0,561.0,300.0,\n",
    "            632.0,540.0,406.2,313.8,576.0,587.6,318.5,400.3,590.6])\n",
    "n=len(x0)\n",
    "lamda=x0[:-1]/x0[1:]  #计算级比\n",
    "rang=[lamda.min(), lamda.max()]  #计算级比的范围\n",
    "theta=[np.exp(-2/(n+1)),np.exp(2/(n+1))] #计算级比容许范围\n",
    "print(rang)\n",
    "print(theta)\n",
    "x1=np.cumsum(x0)  #累加运算\n",
    "z=0.5*(x1[1:]+x1[:-1])\n",
    "B=np.vstack([-x0[1:],-z,np.ones(n-1)]).T\n",
    "u=np.linalg.pinv(B)@np.diff(x0)  #最小二乘法拟合参数\n",
    "print(\"参数u：\",u)\n",
    "sp.var('t');\n",
    "sp.var('x',cls=sp.Function)  #定义符号变量和函数\n",
    "eq = x(t).diff(t,2)+u[0]*x(t).diff(t)+u[1]*x(t)-u[2]\n",
    "s=sp.dsolve(eq,ics={x(0):x0[0],x(5):x1[-1]})  #求微分方程符号解\n",
    "xt=s.args[1]  #提取解的符号表达式\n",
    "print('xt=',xt)\n",
    "fxt=sp.lambdify(t,xt,'numpy')  #转换为匿名函数\n",
    "yuce1=fxt(np.arange(n))  #求预测值\n",
    "print(yuce1)\n",
    "yuce=np.hstack([x0[0],np.diff(yuce1)])  #还原数据\n",
    "plt.plot(x0)\n",
    "plt.plot(yuce)\n",
    "epsilon=x0-yuce[:n]  #计算已知数据预测的残差\n",
    "delta=abs(epsilon/x0)  #计算相对误差\n",
    "print('相对误差：',np.round(delta*100,2))  #显示相对误差\n",
    "'''"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.12101015e-02 2.64159957e+05]\n",
      "[ 1.02314919 -1.10333706]\n",
      "[1.15175675e-01 5.25983987e+02]\n",
      "[ 0.01425587  1.89860155 -4.73846713]\n",
      "2007年人口预测值 284177.1901956999\n",
      "2007年粮食消费预测值： 260900.6772277914\n",
      "2007年畜牧产值： 1114.6587848855643\n",
      "2007年粮食亩产： 1288.9623372343722\n"
     ]
    }
   ],
   "source": [
    "#15.3\n",
    "\n",
    "#灰色G(1,M)模型\n",
    "import numpy as np\n",
    "from scipy.integrate import odeint\n",
    "#导入已知数据\n",
    "x10=np.array([269751,271016,270202,272770,277314,282783])\n",
    "x20=np.array([17675,250172,265415,239894,255345,260253])\n",
    "x30=np.array([566.9,756.0,761.3,681.0,590.0,1236.0])\n",
    "x40=np.array([812.0,1136.0,1108.0,1008.0,1177.0,1251])\n",
    "#求累加序列\n",
    "x11=np.cumsum(x10);x21=np.cumsum(x20)\n",
    "x31=np.cumsum(x30);x41=np.cumsum(x40)\n",
    "z1=(x11[:-1]+x11[1:])/2;z2=(x21[:-1]+x21[1:])/2\n",
    "z3=(x31[:-1]+x31[1:])/2;z4=(x41[:-1]+x41[1:])/2\n",
    "#最小二乘法拟合参数\n",
    "B1=np.c_[z1,np.ones((5,1))];u1=np.linalg.pinv(B1).dot(x10[1:]);print(u1)\n",
    "B2=np.c_[z1,z2];u2=np.linalg.pinv(B2).dot(x20[1:]);print(u2)\n",
    "B3=np.c_[z3,np.ones((5,1))];u3=np.linalg.pinv(B3).dot(x30[1:]);print(u3)\n",
    "B4=np.c_[z1,z3,z4];u4=np.linalg.pinv(B4).dot(x40[1:]);print(u4)\n",
    "def Pfun(x,t):\n",
    "    x1,x2,x3,x4=x\n",
    "    return np.array([u1[0]*x1+u1[1],u2[0]*x1+u2[1]*x2,\n",
    "                     u3[0]*x3+u3[1],u4[0]*x1+u4[1]*x3+u4[2]*x4])\n",
    "t=np.arange(0,7)      #创建合适数量的数组\n",
    "X0=np.array([269751,17675,566.9,812.0])    #导入妹子数据的开始数据\n",
    "s1=odeint(Pfun,X0,t)\n",
    "s2=np.diff(s1,axis=0)\n",
    "xh=np.vstack([X0,s2])\n",
    "pre=xh[-1,:]\n",
    "print(\"2007年人口预测值\",pre[0])\n",
    "print(\"2007年粮食消费预测值：\",pre[1])\n",
    "print(\"2007年畜牧产值：\",pre[2])\n",
    "print(\"2007年粮食亩产：\",pre[3])\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2.87399999999995, 6.04939088432833, 9.36536154002487, 12.8388326140183, 16.5579999999997] \n",
      "----------\n",
      " [2.87399999999995 3.17539088432839 3.31597065569653 3.47347107399341\n",
      " 3.71916738598140] \n",
      "---------------\n",
      " [5.50670620214078e-14 0.102609115671610 0.0210293443034657\n",
      " -0.0834710739934139 -0.0401673859814040] \n",
      "-----------------\n",
      " [1.91604251988197e-14 0.0313023537741337 0.00630187123268376\n",
      " 0.0246227356912725 0.0109180173909769]\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAERCAYAAAB4jRxOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAmUlEQVR4nO3deZzN1f/A8dcbYx97tsRQ2YtEu9CixVeplBYxCaGNipKIX0Khkqi0EEpFlqSVGiHJmqQNY1+yRGKMWd6/P86dxZjl3pk7c+/MfT8fj/u4cz/r+dyZed9zz+ec9xFVxRhjTMFXKNAFMMYYkzcs4BtjTIiwgG+MMSHCAr4xxoQIC/jGGBMiLOCbfEFECgfw3IVERAJ1fmP8xQK+CXoi0hz4WURKBagIozwPY/I1C/gmqIlIDWA+8JGqHgtQMaYDD4jI/QE6vzF+YQHf5JiINBSRr0TkPxE5KCLviki5VOujRERFZGaqZeEikuhZHpnJ4ScC36rqc6n2DRORoSKyTUTiPed8TURK5PA6+ovI1rTLVXU9cCcwzvMB5BciEuHF9RvjN0UCXQCTv4lIBLAE+Bt4EjgTeAyoDlyfZvMmqX4+H8i0XVxEbgAuA+qnWfUq8AAwA4gCWgMPAhWAu32+CHeuO4EXgO3prVfVL0VkHjAW6JSdc6RjP3Av8IOfjmdMpizgm5x6DigBXKOquwBE5B/gRRG5XFWXpdr2bBEp5WmaaZLOsdLqDYxT1QNJC0SkMdALeFBVJ3oWv+X5RtFJRB5W1YPeFt5zM3YwMATYk8Xmg4E/ROQMVd3v7Tky4nkfpuf0OMZ4y5p0TLaJSDHgFuDrpGDv8a3n+cI0uxQCzvP8nGnAF5EyQFvgw3SO8RTweprl6zzrzvCm7KnUBPoCtwPfZLahqm4BVgO3+ngOY4KCBXyTEw2BUsDPaZavxwX21MH6ByCWlEDfBPguk2PXA46q6l+pF6rqelV9QU/P+ncxcAzY4tMVwCGgsarO8XL7paR8aOVIRm34qZYPFpHFIhIjIhNE5FYR2SUi+0XkljTbzxCRvSJyVERWiMh16Zyvg4isFZFjIvKLiHTz3HuJS7rnIs5DIrJRRE6IyJ8iMkhEwvxxzSawLOCbnKjped4LICJFRaQSUNaz7L9U28YBvwBNRCSppr8qk2NXB3Z7UwgRuQC4Cpiqqid9uQBVPaqqXp3HY7enbHnhadyH4o9AH+AV3H2GWFzzEp4b1QuBq4GXgSdwH3zzROTcpAOJyIXALGAn8CiwBngH2AD0BGI8m77kOc7XuPsiX+Ga7V7Ltas0ecYCvsmJ0p7nWM/zTbgbkUmPAWm2X42r2Z8DlPS8zkhxUoJQhjw1z0nAEWCYtwXPgRjcPYu8MEVVhwIfeF73V9VXgcW4G9QA5+KCdqTnm8+buJvKxYBrUx3rSqAwcLeqvq2qXXHfbkqo6mRVjRWRc3AfBi8Aw4F5uPf0E+A+ESmNydfspq3JiUTPc9Lf0fekBJn02sNX43rRXADEc3pTUGp7gWpelGEk0By4U1X3ebF9TlUn65u7/pL0gRiXweukLqMdRKS4iLTEvRdJzT2VUx3rB9zv634RmQVcDpQHUjeZXY3rOTXI80irLu6bgcmnLOCbnEjqPXMGgKr+DSz03MxNz2ogHBeQNgInMjn2VqC6iISr6tH0NhCRO4DHgfGq+pHvxc+WesCveXSu+CxeJ904fwnohvt//h3XDNQyzaa/AL8BY3BNNgAzgQmptqnkee6Be//T8vX+iAky1qRjcmK957lpmuUZtXFvAE7imhwya85BVbfhgtdN6a0XkYuBKbgeQY95VdocEpGSwI24kb/B4glc+/5jQFlVPU9VH0lnuxdxAb8qcAVQS1XvSHPPI6k76zZVXZj0wPWA+g+T71nAN9nmaUL5HvifiNRJter2DLY/iatpQhYB3+MD4MG0ictE5Hzgc9wgqdtV9bSaby7pBuxU1WBq1rgU+E9VX1fV4wAi0j2d7RoDJ1T1gKouU9X0BpgtBBQ3GCy1YcByoIwfy20CwJp0TE71xY20jRKRsUAVXK0zI6tx/fO9Cfiv4gZf9cDdmEVEigCzcTctXwVuTPN58IOnvzwi0gGIUdWvvL+c9IlIVdyNzM45PZafrQfaicibuPf0GqCDZ13xVNt9Dwz0DIr7BfgXd9N2TdJANVXdJCLjgL6ebpqfAw1wA93eyeBDwuQnqmoPe+ToATTD9RyJAbYBj+C6/w31rI8Cojw/98S1RZcAInA1yshMjt0eOApc6nl9gWefjB6RqfbdCvzow3VMAbams7wksAyY6ef3Ld3rT7sciPS8jkhbTs/7+AYutcVxz3t9Je7bz1qgkGe7hrieTLs82yW9XyeBzqnOLcDDuOafE7ibuk8DYYH+O7NHzh/i+SUbE7REpD/uW0Nt9TRbeLnfHcBFqprZNw5vjvMa0AJo48v5g4UnrfR+3M3cebgP5kK4b2N9gaWqekfACmjyjAV8ky+IyPnquiD6ss9M4FH1bWBVesepCiSoH/LnBIqI9AG648ZAlMbV8ncAi4Dhqro3gMUzecQCvjHGhAjrpWOMMSHCAr4xxoSIoO2WWalSJY2IiAh0MYwxJl9ZvXr1AVVNN0140Ab8iIgIVq3KLJmiMcaYtERkW0brrEnHGGNChAV8Y4wJERbwjTEmRFjAN8aYEGEB3xhjQoQFfGOMCRFB2y3TG7GxsRw6dIijR4+SkJAQ6OIYEzSKFi1KpUqVKFu2bKCLYoJIvg34sbGxbN++nfLlyxMREUFYWBhp8qIbE5JUlZiYGHbu3EmxYsUoXrx41juZ4KAK0dFQp07W22ZDvm3SOXToEOXLl6dSpUoULVrUgr0xHiJCyZIlqVSpEvv359sEn6EpJgY6dnRBPxfk24B/9OhRypSxGdeMyUh4eDgnTmQ2T7wJGr//DidOQMmSMH8+1KyZK6fJtwE/ISGBsLCwQBfDmKBVpEgR4uPzarpfk23vvQfNmsGQIe71mWdC4cK5cqp8G/ABa8YxJhP2/xHkjh+Hbt0gMhIuvhgeeyzXT5mvA74xxuRLf/7pgvyUKTB4MCxcCFWr5vpp820vHWOMybdU4dgx+OILuO66PDut1fDzia1btyIibN26NdBFyVLr1q0ZOnRowPbP7ePndvlMARUTA++844J9vXqulp+HwR6shh90oqKi2Lp1K5GRkYEuSra9+eabhIeHB2z/3Bbs5TNB6K+/4Pbb4eefoXFj15xTJP3wu2cP3HknfPSR/1t5rIafxp490KoV7N0bmPNHRUUxZcqUwJzcT+rVq0f16tUDtn9uC/bymSAzcyZceCHs2AGffeaCfSYGDIClS+H//s//RbGAn8Zzz+Xem22MCTHPPAN33AGNGsHatdCuXYabligBIjB9OiQmwuuvu9clSvivOAUu4PftC61b+/4oXNi9ua+/fuqbXbiw78fq29f3ckdGRiIiDBs2jMWLFyMiiMhptf3Dhw/TsWNHSpcuTf369Vm+fHnyuqNHj9KjRw/Kly9P1apVeeyxxzh58qTXZfjnn3/o3Lkz4eHhVK1alWHDhqGqyeunTJlCREQEMTEx9OvXj+rVq7N48eLTjpNRG/eJEye4//77KVeuHOeddx7PP/88tWvXZvz48Vnun3QPY926dRle/65du+jQoQNly5alatWq9OvXj8TERK+v31vpla9169Y88MADtGjRggoVKrBgwQIuvfRSypUrx9y5cwE4cuQIXbt2pWLFilSsWJGuXbty7Nix5GOoKk8++SQVK1akTp06jBgxgoYNG9K/f38A4uLiGDBgAFWqVKFChQpERkby77//+v36jB9deqnrbrl4cZaDqaZOPbWVp2RJuOce/w66LXABP7suuggqV4ZCnnekUCH3OotvX34zdOhQVq5cSY8ePWjWrBkrV65k5cqVtG/f/pTt7r33XqpVq8bcuXMpVaoUffr0SV7Xo0cPvv/+e6ZNm8b48eOZNm0aQ5IGc3jhzjvvZOXKlbz//vsMGzaMMWPGMGrUqFO2UVVuvfVWNmzYwMCBAznnnHO8Pv7o0aOJiopi5syZtG/fnpEjRzJ58uTTrjEzmV3/Pffcw44dO5g9ezYTJkxg6tSpTJ061etj59SMGTMYOHAgNWvWpEOHDjzwwANccsklTJo0CYBHH32UJUuWMG3aNKZPn87SpUsZPXp08v7Tp0/n7bffZsqUKfTt25fBgwczcuRIunfvDsDgwYOZOnVq8u926dKl9OzZM8+uz3hp7lxIqsS0awdjx0LRopnu8sEHLriHh7uKZvHibuBtmTJ+bsdX1aB8XHjhhZqZjRs3Zro+O3r1Ui1USLV4cffcu7ffT5GlZ599Vlu1anXa8ujoaAX09ttvT142a9YsLVKkiKqqbtq0SQGNiopKXj906FCtUaOGV+ddunSpArp27drkZS+99JKWKlVKT5w4oaqqkydPVkA7deqkiYmJGR6rVatW+uyzz562vF27djpgwABVVT1+/LgCunz5cq/2z+r6ExMT9c0339TffvtNVVXj4uK0TZs22qtXL6/L56309m/VqpX27NlTVVW7dOmibdu2VVXVIUOGJP8+p0+frj/99FNyebt06aLXX3998jEefPBBveOOO5JfV65cWWfMmKGq7v0qVqyYTpkyJXn9lClTNCwsTGNiYjIsa278n5gMxMaq9uunCqqXXKIaF5flLomJqi++6HZp1Ur1f/9T7dNHdd0693zLLb4XA1ilGcRV66WTyr590KsX9OwJkya5G7jBJnWNrmLFislD5zdu3Ai4poW0Tp48SdEsahhr166lbNmyNG3aNHlZ69atOXbsGJs2baJRo0aAS7s7duzYbI3ibNCgAUuXLuXw4cMsWrSIIkWKUMfHrIAZXb+IcMcdd/Duu+/y+OOP8+OPP/LPP/9QM5dykqSnWrVqyWVJ/XOS2267jffee49Ro0bxww8/sG/fPq688srk9Q0aNODLL79k7969bNu2jYMHD1KvXj0ANm3aRGxsLJGRkaf14Nq+fTt169bN5aszmdq2DTp1ghUr4OGHYfToDHvhJElIgH793JeBO+5wTTrFiqWsnzDB/8W0gJ/K7NkpP+fGm+0PtWvXTne5etraly1bRsmSJU9ZVySLP7wkaYN40uukY4MLameeeabX5U2tSZMmTJw4kfLlyxMWFsa4ceOoXLmyT8fI6PqPHj1Ks2bNqFy5MnfffTeDBw/m9ddfz1Y5c0NCQgKtW7fm33//5b777uPhhx9m4cKFLF26NHmbJk2asGfPHqpVq4aI8NRTT3HBBRcAKb+DmTNnntaMlpcfaiYd//4LLVq4NpiPP3bdL7Nw4gR07gyffOKa+EePTmlOzk0W8INM8eLFM014VTiDpEpJNfDExMTkWvry5ct5/fXXefvtt7Os4Tdt2pTDhw/z888/06RJEwAWL15MqVKlOPfcc7NxJadSVR566CGWLVtG6dKlqVy5craynWZ0/YsWLSI6OprVq1dTvnx5AB5//HG/lN0ffvnlF1asWMGaNWuSg3hS236SPn368MEHH9CkSRPCw8OpWLFi8rpzzjmHokWLcuLEieTf719//cWwYcMYPXp08jcKk4dUXYN7mTLw4otw+eXgxd/boUPQoQMsWQIvveRq+XnFbtoGmRYtWrBmzRpmz57N999/z6uvvurVfmeffTadOnXi/vvvZ/bs2XzxxRd0796dgwcPZhnsAa644gquvfZa7rzzTubPn89bb73FkCFDGDRoEMVSf8/MJhGhUKFCTJgwgV27drFlyxZ2796d4+MmSQqOkydP5ptvvuGWW27hhx9+CJpskeXLl0dEmDFjBlFRUXTv3p0ZM2acUr7ChQszefJkNm/ezPbt29m2bVvyupIlS9KvXz/69+/PtGnT+Pbbb7n//vv5+eefqZoHOVhMGjt3ugE7X3zhXkdGehXst2+HK65wLT8ffpi3wR4s4Aedq6++mqeeeorevXtzzTXX8NVXX3m971tvvUXLli25//77ueeee2jevDnTpk3zev+PPvqIZs2acdddd/HMM8/Qr18/nnzyyexcRrruu+8+Zs6cSbt27bjgggs488wzadmyJTExMTk+dsuWLRk0aBAjRoygY8eOlC1blgceeIAVK1YExfSXtWrVYvz48UyfPp127dpx4MABBg0axIYNG5K7Vnbr1o1FixZx6623cuGFFxIREUGjRo3Yt28fAM899xz33nsvjz/+OB06dEju/mlZMfPYV1/BBRfAmjXw339e7/bzz3DJJbB7tztEp065WMaMZHQ3N9CPQPTSMbln4cKFWqFCBZ03b57++OOPumzZMn3hhRcUSO5ZE8r+/PNPLVq0qH7wwQf6448/6vLly/Wtt97SsLAw/eKLL7J9XPs/8aP4eNVnnlEVUW3cWNWHv9uFC1XDw1Vr1FD95ZdcLKNaLx0TBFq0aEGHDh148MEH2b9/P4ULF6Zu3bpMnDiR+vXrB7p4AVe7dm0eeughBg0axB5P97DatWszZMgQrrnmmgCXzgCuf/3w4S6H/fjxbmSUF95/H+67z+VL++ILqFEjd4uZGdFUPTCCSfPmzXXVqlUZrv/tt99o0KBBHpbImPzH/k/84NAhqFDB3aSNioI2bbzaTdX1vnnySTcCf84cKFcuNwvqiMhqVW2e3jprwzfGmPQkJLikWnXquGyXIl4H+4QEePRRF+w7dYIvv8ybYJ8Va9Ixxpi0/v7b5TpYuBDuvRd8yI4aE+P62M+eDY8/7nps5kUfe294XQwR6SAiUz0P75OfpOxfXEQ2ichIX/c1xpg88/330LSpS5v71ltukvFSpbza9dAhuPZa13zz8sswZkzwBHvwMuCLSB9gMrAPiAXmiUgXH8/1JFAKGOHjfsYYk3c++ABKl4Yff4Tu3V1Tjhe2bXNjr1audJOXZCdrbm7LsklHRMKBUcBtqvqNZ5kCPQCvUhGKSG3gKeAhVT2a/eIaY0wuOHgQDhxwXWlefhni4twIWi+tWwc33uiac77+2o3JCkbe1PCL4QL1N6mW7fEs99Y44HfctwRjjAkeP/zgmnA6dnSTYZQo4VOwX7gQrrzSzZ2xdGnwBnvwIuCr6gFVTa7Ji0gloDMwz5sTiMj1QHvgADBdRAaLiHcNYsYYk1tUXa76Vq1cvvr33vO5wX36dLjhBoiIgOXL3cRWwcynqxORicAfwHrgRS93S7pJWwWoBjwLLBOR0ybuEpGeIrJKRFbt37/fl6KZbIiKikp3WH56KXiDgYgQFRWV6+cJ1us3fnT0qMtg9sQTcNNNLk1Cs2Ze764KL7zgOvC0bOkSoQVyQJW3fO2WuQw4B2gDXAIsyWxjEWkBNAXeUtWenmVXAt8CXYE3Um+vqpOASeAGXvlYNuMn6U1P6K1169YRFRVF32zescrp/v6Qk+s3+UTx4i6t8bhxLn+9D/mIkvrYT5gAd90Fkyefmsc+mPlUw1fV91W1LfAZMD6r7YGk9HEvpTrG98BfQBNfzm3yTkREBBEREdnad926dbzyyivZPndO9/eHnFy/CWKqbmajAwcgLAwWLYJHHvEp2MfEuHT3EyZA//6uSSe/BHvwIuCLSFEROSvN4gWANwlQklLJpZ2G9wTg/ezaxhiTE0eOuEj9wAPw5ptumY/t9QcPwjXXuJQ6r7wSXAOqvOVNcdsAv4hIhVTLzgW2ZbB9aqsAJVVtXkQqA/WAlT6U0zetW5/+mDjRrTt+PP31U6a49QcOpL/+o4/c+h070l8/f75b/8cf2SpyUnv6xIkTqVy5MlWqVKF///7ExcWdtk1CQgLPPfccERERp6Q/jouLY8CAAVSpUoUKFSoQGRmZnHoXYPfu3dx4440UL16cxo0bs2bNmnTLklEb9vHjx+nTpw8VKlSgUqVKdO7cmQMHDgCuGUREuO+++9i2bRsigoh43Tzi7f6xsbE88MADlClThlq1ajFvXkrfgayu31vpXX9kZCQdO3bk+uuvJzw8nBkzZtC+fXvCw8OZ4JkeLS4ujr59+1K1alXKli3LzTffnJzaOMkrr7xC1apVqV69OsOGDaNFixZ0SpUnd/To0dSsWZMyZcrQoUMH9u7d63P5TRpJ7fNz57qRUE8/7fMhtm51fexXr3aTWj36qN9LmSe8CfjfAjuBBSJyvYhEAk8ArwKISD0RiUhvR1XdDcwAponILSJyI653zz5gVs6LX/CMHTuWd955h5EjR/LGG28wfPjw07bp06cP8+fPp2/fvrRo0SJ5+eDBg5k6dSrjx49n2rRpLF269JQ5YDt37szmzZuZNWsWffr0YfDgwT6VrVevXsybN48JEyYwbdo01q1bR9euXQE31+zKlSt59tlnqVatGitXrmTlypWnnD8z3u7/xBNPcOTIEWbPnk3Dhg25//77SUxM9Or6c2ru3LnceeedXH755dx99920bNmSjh07Jgf8559/nqlTp/Laa68xZ84cdu/efcp8AkuWLOHJJ5/kpZdeYsyYMTz//PP06NGDQYMGAfDGG28wbNgwhg4dmrx/hw4d/Fb+kPTZZ3DppXDypBtB+/jjPjXhgOtjf+mlbs7rr792vTfzrYzyJqd+ALWAT4F/gc24fvlJ66KAKZnsWwzXo2cHcAxYDNTP6pyhlg//u+++U0C//PLL5GVDhgzRihUrnrZNy5YtNTY29pT9jx8/rsWKFdMpU6YkL5syZYqGhYVpTEyM/vrrrwrookWLktf369dP3Z/Aqbp27apdu3Y9ZdmWLVtURHT+/PnJyxYuXKidOnXShISE5GWTJ0/WWrVq+Xz93uwP6MUXX6yJiYmqqrpq1SoFdOfOnVlevy/Su/6uXbtq27ZtVdX9XurWrauqqu+++25yeefNm3fK+zt48GCtX79+8uvRo0frRRddlPz6oosu0pEjRya/Puuss3To0KHJr6OiohTQTZs2+VT+1Ara/4nP/v5btWtX1f37s7X711+rli6tetZZqhs2+LdouYVM8uF71QKlqttU9SZVLaOqZ6vqa6nWtVbVyEz2jVXVAap6lqqWUtVWqvq79x9JoeWyyy5L/rlFixYcPHiQw4cPn7LN2LFjT5u2cNOmTcTGxhIZGZncHBIZGUlcXBzbt2/nr7/+Sj5mkiuvvNLrcm3YsAFVpWXLlsnLrr76aj788EMK5WFDZvfu3ZO7kiZNaxgXF5fl9ftD0ryxInLKz0natWvH7t27ufvuu4mIiGD48OGnzObVoEED/vjjDzZt2sSmTZv4448/qFevHuAmYd+xY0dy05aI0Lp1awA2b97sl/KHjPXrXQL6+Hg44wzXXFupks+HmTbNjZ6tUyd/9LH3hmXLDDKaan6CpKaKtAE1ddBOu9/MmTM555xzTllXs2ZNfv3119OOldGE4N46ceIEq1atSp50Oy/Url073eVZXX9e6NSpEz/99BM9evSga9euREdHM2rUqOT1DRo0QESSJ1bv3LlzcpNNUvnHjRt32gdxnTp18qT8+Z4qvPsuPPQQlC8P0dFezTOb3mFGjXJN/Vdd5bJeli2bC+UNgHx2j7ngW7p0afLPP/30E5UrV6aMF8O8zznnHIoWLcqJEydo2rQpTZs2pVSpUowZM4Z//vmHs88+G+CUG7XLli3zulyNPNWb1OVbvnw5LVu25J9//kleVrx48RxNHJ7V/hl9SGV1/bnt8OHDfPLJJ4wbN47Bgwdz3XXXnfbNom/fvgwdOpRdu3axZ88epk2blvwNoUyZMtSoUYN///03ufwRERGMGTOGrVu35nr5871jx6BrV5fs7IorXMN7NoJ9QoL7vHj6abj7bjdDVUEJ9mA1/KDz6KOPEh8fz/79+3nllVd46qmnvNqvZMmS9OvXj/79+6OqnHnmmQwdOpR//vmHqlWrUq1aNS699FJ69+7NCy+8wM6dOxk/3puhFE6dOnW499576dWrF2PGjCE8PJyBAwdy7bXXnlKDbtasGXv37uWdd96hfv36LFmyxOtryMn+WV1/bitZsiQlSpRgzpw5VKxYkQULFjB27NhTzl24cGFmzpxJ3bp1qVy5MseOHaN27drJ37oGDRrEgAEDqFChAo0aNWLcuHEsX76ciUk9zEzG7rrL3aAdNgwGDXKJbXwUE+OC/Ny5ro/9qFH5r9tlljJq3A/0I1Rv2k6ZMkWrVKmi4eHh+uCDD+rJkydP2yYjJ0+e1P79++sZZ5yh4eHhevPNN+u2bduS1+/atUtvuukmLV26tNatW1eHDx/u9U1bVdVjx45p7969tVy5clqpUiXt2rWrHjhw4LTtJk2apDVq1NAiRYpo48aNfXwnMt4f0O+++y75dXR0tAIaHR3t1fV7K6ObtknLnn32WW3VqpWqnnqTedasWVq7dm0tVqyYtmnTRseOHauFChXSLVu2qKrqp59+qiVKlNDy5curJ+OsnnnmmadM4v7iiy9qjRo1tGTJknrVVVfp+vXrfS5/agXt/+Q08fHuec0aN1N4Nh04oHrZZW5+8nHj/FS2ACGTm7Y2p22QiIqKok2bNgTr78PkzH///ccZZ5zB6NGjufDCCylUqBDR0dH06dOHESNG0KtXr1w5b0H7P0kWE+NSIoSFweuv5+hQ0dEuAdrWrW7C8dtu808RAyWzOW2tSceYPFC6dGmeeeYZXnvtNXbs2EF8fDw1atSgW7dudOni61xCIe6PP9yo2V9+cc03qj73rU+yZo3riXPypEtzfMUVfi5rkLGAHyRat25ttfsCbtCgQcmDrEw2ffAB9Ozpkp998QVcf322D/X11642X6ECfPcdFMQvQmkVtFsSxpiC6u+/XS6cpk1dL5wcBPupU6Fdu5Q+9qEQ7CGfB3yrERuTsQLz/7Fnj2u2qVwZFi921fFsJp9XhREjXA/OVq1ctoXq1f1c3iCWbwN+0aJFTxnFaIw5VUxMDGFhYYEuRs588gnUrw9vv+1eN2vmbtRmQ0ICPPiga/a/5x74/POC1cfeG/k24FeqVImdO3dy6NAh4uLiCk5txpgcUlWOHz/Orl27qFy5cqCLkz2xsS5XfceOrr2lbdscHe74cdde//rr8OSTrkknTXaSkJBvb9qWLVuWYsWKsX//fg4ePJij0Z3GFDRhYWFUqVLFq1HaQSc6Gjp1gpUroV8/NwIqB9H5wAFo3x5WrIDx491I2lCVbwM+uGH4Z52Vdm4WY0y+9ttvsGkTzJnj5p3Ngehod2932zaYNQtuvdU/Rcyv8nXAN8YUEHFxsHQptGnjOsZv2QLlyuXokKHWx94b+bYN3xhTQOzY4brMXHstJKWCzmGw/+ord8hixWDZMgv2SSzgG2MC5/PPXb/6DRtcXgNPVteceO89+N//3KFCqY+9NyzgG2MC45ln3Oins85yk8Wmmts3O1Th+echMtJNMx1qfey9YQHfGBMYpUq5NAnLl2crd31q8fHQu7f7DOncGRYsgPzYQSm32U1bY0ze+eYb93zttfDUU9lOepba8eMuHf6nn7pDjhjhl8MWSFbDN8bkvoQEGDIErrvOtbvkIMNlagcOwNVXw/z58NprMHKkBfvMWA3fGJO79u51U0l9951rYJ8wwS9RecsW18d+xw6XgeGWW3Je1ILOAr4xJvfs3g2NG8OJE26C8fvu88thV692fezj410f+8sv98thCzxr0jHG+Nf+/a5BHVw3md69XZoEPwX7L790fexLlHB97C3Ye88CvjHGP3bvhsceg4gIuPNOOHLELX/+eWjUyC+nmDLF9bE/91zXuad+fb8cNmRYwDfG5Mzu3S7vcJ068OqrLsPlmjV+zT2sCsOHuy8JV13l0uJXq+a3w4cMa8M3xmRPQgIULgz//QeTJ7sbsgMGuMDvR/Hx7vNk0iS4916XGj8UUxv7gwV8Y4xvNmxwnd1jY133mLp1XS0/h/lv0nP8uGsdmj8fBg50rUPW7TL7rEnHGOOd1atdfuHzznMR+NxzITHRrcuFYL9/v2u+WbAAJk60AVX+YDV8Y0zW3nvPNdmUKwfPPutmo6pQIddOt3kz3HBDSh/7HKbFNx4W8I0xp1OFb791+YWvuMIlORs5Evr0yfUkNatWudPFx8OiRXDZZbl6upBiTTrGmBSqrg3lssvgmmvghRfc8kqVXKKaXA72X3zhMl2WLAk//GDB3t8s4BtjnK++gmbNXEf3PXvcjN8zZ+bZ6SdPdnPP1q3r+tjXq5dnpw4Z1qRjTCiLj3e1+rAwl5zm2DEXee+5xy3LA0l97IcMgbZt3dyz4eF5cuqQYzV8Y0LRyZPwzjtuqOrbb7tl3bu7CcQjI/Ms2MfHwwMPuGDfpQt89pkF+9zkdcAXkQ4iMtXzaJ+dk4lIUxGJE5GI7OxvjMmhmBiXR/icc1yAL1cuZVrBsDA3kCqPHDvmMly+9RYMGuTSJuTR50zI8irgi0gfYDKwD4gF5olIF19OJCKFgElYM5IxgXP77fDww1CzprtDunKla0fJY0l97D//3N0qGD7c+tjnhSyDr4iEA6OA21T1G88yBXoAU30414OApToyJi8dPuxGLfXs6XraDBwI/fvDlVcGLMJu3uzy2O/cCbNnw803B6QYIcmb2nYx4KGkYO+xB2jq7UlE5EzgeWAg8JovBTTGZMOBA/DKKzB+PPz7L9So4RrJA5xLeOVK18c+MdF187/00oAWJ+Rk2aSjqgdUNbkmLyKVgM7APB/OMx74HFjgcwmNMd5LTHQ1+Fq1XC6Ctm1d5souPrXA+tWePS5//fTpro996dIuj70F+7znU3u6iEwEOgHfAy96uc9NQCugAVAyi217Aj0Batas6UvRjAlt//wD5ctDoUIQHe1y3gwcCA0bBrpkPPccLFkC33/vuvkvWABVqwa6VKFJVNX7jUXuAboCFwHtVXVJFtuXBjYCz6rqZE/vnGigtqpuzWzf5s2b66pVq7wumzEhadMmGDUK3n8ffv7ZjVpKTHSBP4Di411NPjb29HXFi7vOQiZ3iMhqVW2e3jqfaviq+j7wvohMxzXTNM1il+HAn6o62ZfzGGOysHGja7KZMcP1ZezRI6UDex4G+9hY+PNPV5zffnPPGze6ZXFxp25booT74jFmTJ4Vz6ThTS+dokAVVd2RavECoKMXx+8A1PL06kktWkTeU9VIbwtqjPE4fBiaN3eB/fHH3bSCudxGcuwY/P776YF98+aUDMkirkt/gwYuO0PDhjBvHsyd6yYsiY11qXisOSdwvKnhtwE+EpE6qnrIs+xcYJsX+94IpJ6bpjruw6IdsMGXghoT0pYvdznoR4xwg6U+/tjd9axY0a+nOXzYBfTUQf2332Dr1pRtihRxqfDPP99NTtKwoQvydeu6Wnxqc+dCr16uV+ikSe4GrgmcLNvwRSQMWAscBYYBVYFXgYGqOkFE6gGxWbXJe44VgbXhG+MdVTd56/DhLk9wxYqwfj1Ur57jQ+/ff2pQTwrsu3enbFO8uEtg1rBhSlBv2NAN0rURscErR234qhonIu1wbfYfA/uBp1V1gmeTN4GtQKRfSmuMcYnMunRx/RerVoWxY13SmVKlvD6EqgvgaWvrGze6bvpJSpVygfzaa1OCesOGEBGRp5kWTB7w6qatqm4DbspgXWtvT+ap1dsAamPSk5gIu3bBWWdBlSpw4gRMmADdurnqdia7bd9+elDfuNGNuUpSvrwL5B06pAT1hg3dmCxLaxAaLK+NMYEWH+/a5EeMcFksf/vNVbtXrjwlEsfHu4p/2sD+++9usu8kVaq4QN6586nNMVWqWGAPdRbwjQmUkyfd8NORI11/+oYNYehQYmPhry2wcaOcEtj//NPtkuSss1wg79kzJag3aOD3+7imALGAb0yAxM6aT7H77+dgRDO+vG02n8TfzK/DCrG5CyQkuG1EoHZtF9BvvDElsNevn+uzDZoCyAK+MbnsyBFXS/9r3THKffQm2/eXYOzx3myP7kAbvmbh1msoslM45xxo3BjuuCPl5mm9eqd3dTQmuyzgG+OjPXtc//OPPjp1ENGBA6ffNP3tNzi66wgPMoF+vMwZHOCrsrdz8Q29iYwsTMOG1zLO09WxaNGMz2mMP1jAN8YHqvDUUy4ZWMeOcN55KQF+//6U7UqVcrX0QRHvc9/BByl+4gjHWt1IwnODuK7lZVwXuEswIcwCvjEZOHHCBfOff3bjncaNS0kjAK6L/LJlrp29W7eUHjGNKu7lzBpCoWpVYElNeOVqGDSIUs2aBe5ijMECvjHJA5TWr3ePpAD/++8pN0+LF3ft6//9Bzt2uMRgJUq4OVnHjvU07ezYAaNHu0lau3VzfehbtnQPY4KABXwTUk6ccM0vSUE96fngwZRtatZ0eWI6dHDP55/vcscULgy9e7ucMMWLu2RgZctC1WOboecLbhZuVTdCtm/fAF2hMRmzgG8KJFU3aDV1UP/5Z9eXPanWXqKEa4O/5RZo0sQF9vPOcyNSM7JvXzrJwJ5/Hj74wKUoHjDAzTZlTBDyaQKUvGTJ04y3YmLg119PDe7r18OhQynb1KrlAnpSYG/SxKXyzVaumHXr3KjYAQNcmuJdu1xDvh+SmhmTU36bAMWYQFKFnTtPDepJtfakm6klS7pa+m23nVprL1cuhydPSHDpiceNg6goN+rplltcwD/zzBwe3Ji8YQHfBKXjx12tPXVwX7/eTd2aJCLCBfXbb08J7mefnQsTPqnCxRfD6tWugf/FF13zTY4/RYzJWxbwTUCpus4taW+i/vVXSq29VClXS7/jjpTmmMaN3Q3TXPP7765dfuhQ9wnSqxdUqAA33eRmADEmH7K/XJNnjh07tdae9HzkSMo2deqkzKSUFNxr186jaVoTE+HLL+HVV+Grr9zQ19tvd5823bvnQQGMyV0W8I3fqcK2baf3kNm0ya0DKF3aBfS77kppjmncOIAJwTZvhhtucF8tqlWD555zXXEqVw5QgYzxPwv4JkeOHYNffjl90FLqiTfOPtsF9XvuSam1R0TkUa09M5s3u0fbtq4bT6NGrgmnY0dLbGMKJAv4JsNkYKmpuoms09baN29OqbWHh7uA3rlzyoCl885ztfmgoermh331VfjsM5dUPjratcvPmRPo0hmTqyzgG557DpYuhf/7P5g40aUPSKq1p+4lc/So217EZXds0gTuvTelSSYiIshnVPrqK3jsMTfU9owz4Jln3M3YgH/VMCZvWMAPYSVKuFQDSV5/3T1SK1PGBfQuXVKaYxo1CrJae2a2bXPNM9WquVFWxYq5FAidOmU6T6wxBZEF/BC2ZYvLF/PTT+61iGvhuPNOuPxyF9xr1gzyWnt6VN1XlnHjXDPNI4/Ayy/D1Ve7vvT57oKM8Q8L+CEsOtrNkw2usnvyJLRrBy+8ENhy5cgHH8CYMbB2rUuK88QT8OCDbp0FehPiLOCHqN27XfqBEiVcjf6RR1IlA8tv9u93bfLg2ulPnoQ333R3j0uWDGzZjAkiFvBDUGws3Hqruwm7YoXr/w4ufXu+smKFa7aZOdP93KyZu4hSpaw2b0w6LOCHGFXo08fFx08+SQn2+UZcHMya5QL9ihXurvJDD6UMkMo3d5ONyXsW8EPM66/Du+/C4MGulp9vJCS4XjYnTriulFWrwvjx0LWrGwBgjMmSBfwQ8v338Oij8L//uQGl+cLatW6Q1C+/uDvM4eGuW9G551r/eWN8ZP8xIWL7dpcx4OyzYfr0II+V8fGu2ebKK127/MyZLj1xTIxbX69ekF+AMcHJavghICbGzdURGwvz5uVyWmF/+OQT13UoIsLNEN6tm+WeN8YPLOAXcKou6ePatfDpp65yHHR+/dU12zRq5PqH3nKL+2Rq1y6bcxAaY9Jj34sLuFdecU04//d/ru0+aCQmuikDr7nGdRWaOtX1pweXCuGmmyzYG+NnVsMvwBYudANNb7sNBg0KdGnS6N4dJk9288GOGOGmDKxUKdClMqZAsxp+ARUd7fKDNWzocoUFfBzSn3+65prt293r7t1dPuboaBg40IK9MXnAavgF0LFjLilaYiLMnRvAsUiq8PXXbpDUF19AWJjLylazJlx2WYAKZUzo8rqGLyIdRGSq59Heh/2qiMgcEflPRGJEZIaIlMpecU1WVOG++2DDBvjwQ9cNMyDi4qBpU7j+elizxnX8377dfe0wxgSEVwFfRPoAk4F9QCwwT0S6eHmOWUATYCAwEujoeTa54IUXXLf1UaPguuvy+ORbt6Yk1A8Lc71tpk1zgf7ZZzOeTssYkydEk+any2gDkXBgF3Cbqn7jWTYJaKCqLbPY92pcwG+oqns8yyYC/1PVmpnt27x5c121apXXF2Jcq0m7dq4L+/vv51G7vSosXuyabT791A2I2rrV3Yw1xuQ5EVmtqs3TW+dNDb8Y8FBSsPfY41melVXApUnB3uMgEObFvsYHf/0Fd93lJi15++08Cva//OKabdq0gSVL4Kmn3E1YC/bGBKUsb9qq6gFgatJrEakEdAbe9WLfI8CRNIuvBX70rZgmM0ePws03u1aUuXNzOQX8zp2wbx9ceKGbHqtUKfcJc/fdLrm+MSZo+dRLx9Mc0wn4HnjR15OJSFvgYuCqDNb3BHoC1KyZaYuP8UhMdPPN/vknfPMN1KqVCydRheXLXbPNJ5/ABRe4RGblysEPP+TCCY0xucHXfvjLgNVAG+ASX3YUkRLABOALVf0uvW1UdZKqNlfV5mckzWBkMvXcc65W/9JLrmXF7z7/HFq0cN0pv/oK+vZ1d4WNMfmOTzV8VX0feF9EpgPjgaY+7P4iUBH3YWH8YN4819uxa1d4+GE/HnjvXtdUEx7ufj52DCZOhHvvtQlGjMnHsqzhi0hRETkrzeIFQH1vTyIitwIPAT1UdadvRTTp2bjRTdnaogW88YafbtKuWuWCes2arl0eXHvRr79C794W7I3J57xp0mkD/CIiFVItOxfY5s0JRKQp8B7wmqp+4nMJzWkOH3YjaUuVgtmzoXjxHBxM1R3k8svdp8fcuW5GqaRMa0WKWO55YwoIb5p0vgV2AgtEZBhQFXgCN5AKEakHxKrq1rQ7ikgY8BHwDzBDRFL3DV2vqidzVvzQk5DgOsRs3QrffQc1auTwgHFxLuf8vn0utWZkZD5ImG+MyQ5vumXGiUg7XJv9x8B+4GlVneDZ5E1gKxCZzu7nAXU9Py9Ls662Zz/jgyFD3ACrN95wlfJsi49388OWLg1ffun6clo6YmMKNK9u2qrqNuCmDNa1zmS/NUCg8zQWGDNnukzCPXvCAw/k4EAnTrhRWocPu76cNgm4MSHBGmfzifXrXWvLZZe5yaGy7d9/4cYbXVv9rbe6NnpjTEiw//Z84OBBd5O2XDk3t3cxb5JapGf/frjhBli3zk2Ddc89/iukMSboWcAPcvHxLhnarl3w/fdQrVoODnbXXa6LZdJ8scaYkGIBP8g99ZSbqvDdd+Hii3N4sFdfdV8XWmaa5NQYU0BZG34Qe/9912Py4YfdpCbZsmYNDB7s+ts3bGjB3pgQZgE/SK1Z46Z9bdXKBf1sWbwYWreGqVNdzd4YE9Is4Aehv/92N2krV3ZdMcOyM3vAp5+6Ka9q1IBly2yScGOMBfxgExcHt9/uOtTMmQPZSho6bZrrcnn++e5Ob46H4xpjCgIL+EHmscdcjH77bWjWLJsHCQ+Ha66BRYusZm+MSWYBP4i8+y689ho8/ng2usirwtq17ucOHVz+BRtBa4xJxQJ+kFixwmUgvvZaGDXKx50TE11XnhYt3N1eyKNJbY0x+Yn1ww8Ce/a4JvcaNeDDD33MdhAX52ZAmTEDnnjCTT9ojDHpsIAfYLGxcNttLo/Zjz9ChQpZ7pLi+HHo2NE134waBU8+mVvFNMYUABbwA0gVHnrIzQ8+cyacd56PB5gxw6U2njQJevTIlTIaYwoOC/gB9OabrjfO00+7irrXVF0bfbdurgkn2915jDGhxG7aBsiSJe4+6403wv/9nw87RkfDpZe6SW1FLNgbY7xmNfwA2LnT1ehr13b5cryeaOrXX6FtW4iJgaNHc7WMxpiCxwJ+HjtxAm65xcXsqCiX494rK1a4rwPFirmRWY0b52IpjTEFkQX8PKQKvXrBqlUuJX2DBl7uuHo1XH01VK3qpiSsXTtXy2mMKZisDT8PvfoqvPceDBsGN6U7Q3AGGjWCLl1g6VIL9saYbLOAn0e+/dalTOjQAZ55xsudPv4Y/vkHiheHiRNdDd8YY7LJAn4e2LoV7rgD6tVzqekLefOuv/ACdOrkno0xxg+sDT+XHT/uavXx8TB3rhf5zFTdvIYvvujmoPWpz6YxxmTMAn4uUoX774f16+Hzz+Hcc7PYISHB3dV9+23o0wfGj/fy64AxxmTNokkuGjPGJUMbORKuv96LHQ4ehK+/do38r71mwd4Y41dWw88lX33lWmbuuAMGDMhi42PH3I3ZypXh55996JxvjDHesypkLti0Ce68042NevfdLFLTHzwIV13l8iyABXtjTK6xgO9nR4+6m7SFCrmbtKVKZbLxrl1w5ZWuVu9Vm48xxmSfNen4UWKim4vkt99cU3ymY6T++stNb3XokEtx3Lp1XhXTGBOiLOD70YgRMGcOvPSSy4SQoZMn4brrXNv9d9/BhRfmWRmNMaHLAr6fzJ8PQ4ZA587Qt28WGxct6pLhn3UW1K+fF8UzxhgL+P7w++8u0Ddr5iafyvAm7YIFbgLb7t1dc44xxuQhu2mbQ0eOuJu0xYq55pwSJTLY8P333YZvveWG3RpjTB6zgJ8DiYlwzz2weTPMmuVaaNL12mvuK0DLli69cRH7YmWMyXteB3wR6SAiUz2P9j7sV0xEJorIIRH5U0RuyF5Rg8+zz7pWmnHjXO/KdD33nOtjf/PNLr9CmTJ5WkZjjEniVcAXkT7AZGAfEAvME5EuXp7jVeAu4GHgFWCWiNTzvajB5ZNPYPhwlyund+9MNgwLg8hI9xWgePG8Kp4xxpxGVDXzDUTCgV3Abar6jWfZJKCBqrbMYt8awDagi6q+71n2NpCoqj0z27d58+a6atUqry8kL23YAJdcAued56YpLFYszQZxcW64bYMGLoMaZDHc1hhj/ENEVqtq8/TWeVPDLwY8lBTsPfZ4lmelFZAAzE617DMgs17qQe3QIdc6U6aMq+WfFuxjYuC22+DSS2H/fhfoLdgbY4JAlncPVfUAMDXptYhUAjoD73px/OrAFlWNSbVsB1BLRAqrakLqjUWkJ9AToGbNml4cPm/Fx7sU9Tt3wuLFUL16mg2OHHFzFy5ZAhMmwBlnBKScxhiTHp966YjIROAPYD3wohe7FAcOp1kWAxQGyqXdWFUnqWpzVW1+RhAGy6efdikTJkxwTTqn+PtvaNMGfvgBPvggi4Z9Y4zJe752y1wGrAbaAGlDXnpicU06qZ30PGfUYz0ozZgBo0e7eUm6d09ngxdfdCOw5s93qTKNMSbI+BTwVfV9VW2La4cf78Uuf+OadVIr73k+7su5A2ntWtcbp2VLePnlNCuTbsqOGAHLl1vWS2NM0Moy4ItIURFJO6RoAeBNEpj1QE0RSR30LwBOAP94XcoA2r/fDZCtWBFmznRpcJKtXAmtWsGBA25FkyaBKqYxxmTJmxp+G+AXEamQatm5uO6WWVkH7AT6AYhIEaAH8K1m1R80CMTFQadOsG+fS5tQpUqqlYsWuYlLdu50SfCNMSbIeTPG/1tc0F4gIsOAqsATwEAAzyCqWFXdmnZHVU0Ukf7ADBGpC1QDmgGZ9t8PFv37u+zFU6dC89S9WufMce30deu6uQxP665jjDHBJ8savqrGAe2A/cDHwGDgaVWd4NnkTWBoJvt/DLTH9co5DrRV1R9yVOo88N57LmVCv35w772pVsyeDR07utSY6fbNNMaY4JTlSNtACeRI259+crlxLr/cVeBPyXW2d69LfP/yy1nMX2iMMXkvpyNtQ8revXDrrVCtGnz0kSfYq8K0aW7kVdWqLum9BXtjTD5jeXpTOXnStdYcOuR6WFaqBCQkwIMPuhmqkiatNcaYfMgCfiqPPgrLlsGHH3p6WJ486RrwP/4YBg6ELt4mCDXGmOBjAd9j0iR44w148knXFZNjx1zbztdfuyG2TzwR6CIaY0yOWMDH1eofesgNkn3+ec/CzZthxQp45x3o1i2g5TPGGH8I+YC/a5fLZlyrlst5VjjmPyhdGs4/H7ZsgQoVsj6IMcbkAyHdS+fECddqc+wYzJ0L5Q9tdoF+gmeIgQV7Y0wBErI1fFWXwfinn9xYqkYJ6+GK61w+hYsuCnTxjDHG70K2hj9hAkyZ4sZQ3VLlB5cErXBhN3lJixaBLp4xxvhdSNbwo6Kgb19o3x6efWAv1G3rUiR8841rzDfGmAIo5AL+tm1w++1w7rkwfToUKuMZOXvNNVC5cqCLZ4wxuSakmnSOH4dbbnHjqaLuepMyq79zK+6+24K9MabAC5mArwo9esC6tcpPNz9PlWd7wVtvBbpYxhiTZ0KmSeell2DGB4ksv/wJ6k17GTp3hnffDXSxjDEmz4REDf+bb2Bg/3gW1bqfi5e9DI884hLeh4UFumjGGJNnCnwNf8sWlxunfsNCXHF+InQbBoMHg0igi2aMMXmqQAf8//6De9r/y5mJh5kzryZhdaZYoDfGhKwCG/BV4dG79zN+4w3Ur3mc0rXWgxTYyzXGmCwVyAi4Zw/cdtEO3t15LWcX2UbYxFlp5ik0xpjQUyCj4JBOf/DhzmspX+gIRRZ9DVe2DHSRjDEm4ApUwC9RwmXA/IzHKEYsLRMX83OrphQvDjExgS6dMcYEVoHqlrllC9x0E9wn73E5y/irZFPuuQeiowNdMmOMCbwCVcOvVs3lQDsolTharBInT0CZMlC1aqBLZowxgVegavgA+/ZBr17w44/uee/eQJfIGGOCQ4Gq4YObzCRJ0sRVxhhjCmAN3xhjTPos4BtjTIiwgG+MMSHCAr4xxoQIC/jGGBMiLOAbY0yIEFUNdBnSJSL7gW05OEQl4ICfihMK7P3yjb1fvrH3yzc5eb9qqeoZ6a0I2oCfUyKySlWbB7oc+YW9X76x98s39n75JrfeL2vSMcaYEGEB3xhjQkRBDviTAl2AfMbeL9/Y++Ube798kyvvV4FtwzfGGHOqglzDN8YYk4oFfGOMCREFMuCLyH0iEhXocgQ7EakiInNE5D8RiRGRGSJSKtDlCmYi0kFEpnoe7QNdnvxCRJqKSJyIRAS6LMFKRO4UEU3zWOjPcxS4gC8iFwOWCd87s4AmwEBgJNDR82zSISJ9gMnAPiAWmCciXQJbquAnIoVwNyEL3PwbftYU+BRokerRx58nKFC/ABG5CpgD/B7osgQ7EbkaaAw0VNU9nmVVgQ7AIwEsWlASkXBgFHCbqn7jWaZAD2BqIMuWDzwI1A90IfKBpsC3qroqt05Q0Gr4LYEuuE9Jk7lVwKVJwd7jIBAWoPIEu2LAQ0nB3mOPZ7nJgIicCTyP+xZpMncBsDo3T1DQAv5zqjov0IXID1T1iKqm/SZ0LfBjIMoT7FT1gKom1+RFpBLQGbC/t8yNBz4HFgS6IMFMRKoDlYHnPffUdovI8yLi11aYAtWko6qJgS5DfiUibYGLgasCXZZgJyITgU7A98CLAS5O0BKRm4BWQAOgZICLE+wuBBKB74BBuHtrI4HjuG9IflGgAr7JHhEpgbvR/YWqfhfo8uQDy4BzgDbAJcCSwBYn+IhIaeA14AlV/dt652RpCdBEVTd4Xi8SkbJAT/wY8Atak47JnheBirg/LpMFVX1fVdsCn+GaLMzphgN/qurkQBckP1DVw6mCfZLlQE0RKeOv81jAD3EicivwENBDVXcGujzBSkSKishZaRYvwHqfZKQDcHVSf3Ig2rM8WkSmBKxUQUpEzhaRumkWl/c8F/fXeaxJJ4SJSFPgPeA1Vf0kwMUJdm2Aj0Skjqoe8iw7l5xN0lOQ3QgUTfW6Ou4Dsh2QtiZr4FHgLOCWVMvuBrar6t/+OokF/BAlImHAR8A/wAwRST3ZwnpVPRmYkgWtb4GdwAIRGQZUBZ7AuhumS1U3pn4tIoc9P25U1e15X6KgNxlYISKv4rpMtwf+hxvn4TcW8EPXeUDSV8hladbVBrbmaWmCnKrGiUg7XJv9x8B+4GlVtVHdJsdUda2I3IXrmdMT2Ai0U9XP/XkeS49sjDEhwm7aGmNMiLCAb4wxIcICvjHGhAgL+MYYEyIs4BtjTIiwgG+MMSHCAr4xxoSI/welS7/LDO/7SgAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#15.4\n",
    "\n",
    "import numpy as np\n",
    "from sympy import Function,diff,dsolve,symbols,solve,exp\n",
    "import matplotlib.pyplot as plt\n",
    "x0=np.array([2.874,3.278,3.337,3.390,3.679])     #导入数据\n",
    "n=len(x0)\n",
    "x1=np.cumsum(x0)         #计算第一次累减序列\n",
    "ax0=np.diff(x0)          #计算一次累减序列\n",
    "z=0.5*(x1[1:]+x1[:-1])       #计算均值生成序列\n",
    "B=np.c_[-x0[1:],-z,np.ones((n-1,1))]\n",
    "u=np.linalg.pinv(B).dot(ax0)\n",
    "p=np.r_[1,u[:-1]]        #构造特征多项式\n",
    "r=np.roots(p)       #就特征根\n",
    "xts=u[2]/u[1]        #常微分方程的特解\n",
    "c1,c2,t=symbols('c1,c2,t')\n",
    "eq1=c1+c2+xts-2.874    #这里的2.874是输入数据的首个数据\n",
    "eq2=c1*np.exp(4*r[0])+c2*np.exp(4*r[1])+xts-16.558      #这里面的数字4是导入数据减一，而后面的16.558是所有数据之和\n",
    "c=solve([eq1,eq2],[c1,c2])\n",
    "s=c[c1]*exp(r[0]*t)+c[c2]*exp(r[1]*t)+xts         #微分方程的符号解\n",
    "xt1=[]\n",
    "for i in range(5):xt1.append(s.subs({t:i}))\n",
    "xh0=np.r_[xt1[0],np.diff(xt1)]\n",
    "cha=x0-xh0       #计算残差\n",
    "delta=np.abs(cha)/x0      #计算相对误差\n",
    "print(xt1,'\\n----------\\n',xh0,'\\n---------------\\n',\n",
    "      cha,'\\n-----------------\\n',delta)\n",
    "#画出该模型的图像\n",
    "hgq=np.arange(1,6,1)\n",
    "plt.plot(hgq,x0,'b*-',label='the original image');plt.legend()\n",
    "plt.plot(hgq,xh0,'r--',label='predict the image');plt.legend()\n",
    "plt.title(\"GM(2,1) image\")\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}